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ABSTRACT 


The  use  of  computational  analysis  in  the  design  of  propulsive 
nozzle  installations  has  recently  expanded  as  advanced  digital  com¬ 
puters  have  been  developed  which  result  in  lowering  computational 
costs  versus  actual  wind  tunnel  test  costs.  Although  a  range  of 
numerical  techniques  has  been  applied  in  this  area,  only  those  utiliz¬ 
ing  the  full  Navier-Stokes  equations  across  the  flow  domain  have 
successfully  simulated  the  viscous  phenomena  associated  with  aft-end 
flowfields.  Navier-Stokes  methods  are  particularly  useful  for  predicting 
off-design  nozzle  characteristics  where  the  overexpanded  or  underexpanded 
flowfield  is  more  complex  and  where  viscous  regions  are  more  prevalent 
than  at  on-design  conditions.  One  feature  typical  of  these  off-design 
conditions  is  the  appearance  of  a  strong  normal  shock  wave  referred  to 
as  a  Mach  disc.  Viscous  nozzle  flowfields  containing  this  phenomenon 
have  not  been  adequately  simulated  in  the  past.  This  research  details 
the  development  of  a  numerical  Navier-Stokes  method  capable  of  accurately 
predicting  supersonic  coflowing  nozzle  flowfields  which  contain  both 
highly  viscous  regions  and  complex  shock  structures  typified  by  the  Mach 
disc  formation. 

Numerical  solutions  to  the  Navier-Stokes  equations  are  obtained  for 


a  domain  containing  an  axisymmetr ic  coflowing  nozzle  with  a  thick  base 
annulus  (M  =  1.94,  M.  =  3.0,  Re  =  2.2x10^).  Five  nozzle  pressure 
ratio  conditions  ranging  from  a  highly  overexpanded  case  (P^ /Pa  =  0.15) 
which  exhibits  a  Mach  disc  shock  formation,  to  a  slightly  underexpanded 


case  (P./P 
J  “ 


1.59)  are  examined  and  solved  numerically.  The  weak  con- 


xv 


servative  form  of  the  two-dimensional  (axisymmetric) ,  time  dependent 
Navier-Stokes  equations  is  solved  using  MacCormack's  explicit  finite 
difference  method.  This  algorithm  is  an  efficient  Lax-VJendroff  type 
differencing  scheme  of  second  order  accuracy  which  utilizes  time¬ 
splitting  and  two-step  predictor-corrector  techniques.  An  adaptive 
grid  scheme  is  utilized  in  the  wake  of  the  nozzle  base  annulus  that 
allows  the  fine  mesh  region  of  the  computational  grid  to  remain  in 
the  mixing  layer  containing  high  flow  gradients  as  each  solution 
progresses  towards  convergence.  Appropriate  numerical  boundary 
conditions  are  applied  that  allow  the  computational  domain  to  be  re¬ 
stricted  to  a  compact  region  surrounding  the  nozzle.  Locally  depend¬ 
ent  eddy  viscosity  modelling  is  applied  in  the  form  of  a  Cebeci-Smith 
two  layer  model  in  the  boundary  layer  regions  on  the  nozzle  walls, 
and  a  form  of  the  Prandtl  mixing  length  model  in  the  nozzle  wake  region. 

The  numerical  solutions  successfully  reproduced  all  of  the 
essential  nozzle  flow  features  including  boundary  layers,  corner 
expansions,  recompression  shocks,  the  separated  recirculation  region 
along  the  nozzle  base  wall,  and  the  evolution  of  the  near  wake  to  a 
far  wake  type  of  flow.  Correct  transition  from  regularly  reflected 
shock  waves  at  the  line  of  symmetry  in  the  jet  core  to  the  strong  Mach 
disc  shock  reflection  was  numerically  achieved,  as  was  the  simulation 
of  the  subsonic  embedded  region  immediately  behind  the  Mach  disc  shock 
structure.  Mumerically  obtained  nozzle  base  pressure  coefficients  were 
within  seven  percent  of  the  experimentally  determined  values  for  all 
cases  where  the  flow  obeyed  the  assumption  of  remaining  attached  in  the 
divergent  portion  of  the  convergent-divergent  nozzle. 

Present  solution  times  for  2,500  point  grids  are  on  the  order  of 


two  to  four  hours  when  run  on  a  Cyber  750  computer.  A  fully  vectorized 
version  of  the  present  computer  code  can  be  expected  to  converge  within 
five  minutes  on  a  CRAY-1  computer  for  similar  grids,  allowing  the  compu¬ 
tation  of  more  complex  nozzle  geometries  and  better  resolution  in  the 
boundary  layers  through  the  use  of  a  finer  mesh  in  future  efforts. 


CHAPTER  I 


INTRODUCTION 


1 . 1  BACKGROUND 

The  increased  importance  of  the  aft-end  drag  problem  associated 
with  nozzle  installations  in  current  and  future  high  performance  air¬ 
craft  has  led  to  extensive  and  very  costly  experimental  nozzle  test 
programs.  Any  technique  which  can  reduce  this  requirement  for  wind 
tunnel  testing  in  the  design  of  nozzle  installations  will  result  in  a 
significant  savings  to  the  technical  community  of  both  time  and  resources. 

Computational  aerodynamics  shows  great  promise  as  a  field  which 
can  have  a  favorable  impact  on  this  requirement  for  nozzle  design 
information.  Current  Navier-Stokes  techniques  in  this  area  utilize 
advanced  digital  computers  to  simulate  the  flowfield  surrounding  the 
nozzle  at  projected  flight  conditions.  It  has  been  shown  that  boundary 
layer  and  shear  layer  growth,  areas  of  separated  flow,  shock  wave 
formation  and  interactions,  and  jet  plume  blockage  and  entrainment 
characteristic  of  nozzle  flows  can  be  analyzed  using  computational 
techniques.  Unlike  experimental  testing,  computational  analysis  is 
not  necessarily  restricted  by  wind  tunnel  Reynolds  number  or  nozzle 
exhaust  temperature  limitations.  Flowfields  analyzed  computationally 
can  also  eliminate  the  undesirable  effects  of  support  stings  and  test 
section  walls  that  occur  during  experimental  testing.  As  more  ad¬ 
vanced  computers  are  developed,  the  cost  of  numerical  analysis 
decreases.  Since  the  cost  of  wind  tunnel  testing  is  steadily  in¬ 
creasing,  computational  analysis  is  being  utilized  more  extensively. 

Several  of  the  first  computational  solutions  to  include  viscous 
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effects  inherent  to  aft  end  or  nozzle  flowfields  consisted  of  patching 
techniques  that  divided  the  field  into  predominantly  inviscid  and 
viscous  regions.  Grossman  and  Melnik  (1),  and  Cosner  and  Bower  (2) 
obtained  transonic  boattail  nozzle  solutions  using  iterative  techniques 
that  divided  the  flowfield  into  an  inviscid  freestream,  an  inviscid 
jet,  and  a  viscous  boundary  layer  and  mixing  layer  region.  The  free¬ 
stream  solution  procedure  assumed  irrotational  potential  flow  that 
could  be  solved  by  a  relaxation  algorithm  applied  to  the  potential  flow 
equations.  The  rotational  inviscid  supersonic  jet  was  solved  using  a 
hyperbolic  marching  technique.  Imbedded  shocks  in  the  jet  were  ex¬ 
plicitly  fitted  to  satisfy  the  Rankine-Hugoniot  equations.  The  viscous 
mixing  region  was  assumed  to  be  isobaric  and  was  solved  using  integral 
techniques.  Each  region  was  solved  separately  and  patched  together 
iteratively  using  pressure  and  flow  direction  conditions  at  the  common 
boundaries.  Separation  regions  could  not  be  accounted  for,  so  equiva¬ 
lent  fitted  body  blending  was  used  to  obtain  reasonable  flow  solutions. 

Pergament,  Dash,  and  Wilmoth  (3)  introduced  a  displacement  thickness 
correction  to  the  inviscid  plume  boundary  to  account  for  the  effects  of 
jet  entrainment  on  the  inviscid  external  flow  calculation.  Their  analysis 
also  included  the  effects  of  species  mixing  and  pressure  gradients  in  the 
mixing  region,  but  still  could  only  account  for  separation  regions  by 
body  blending  techniques.  Yeager  (4)  attempted  to  include  a  fourth 
separation  region  involving  recirculating  flow  that  was  defined  by  a 
dividing  streamline  which  connected  separation  and  reattachment  points. 

The  extent  of  this  region  was  determined  using  local  control  volume 
analyses,  and  it  was  found  that  reasonable  reattachment  points  could  only 
be  predicted  through  the  application  of  empirical  corrections  during  the 


solution  procedure.  Although  these  iteratively  patched  solutions  gave 
reasonable  results  for  specific  data  sets,  the  required  amount  of 
empirical  matching  and  explicit  fitting  limits  the  use  of  this  type 
of  computational  method  as  a  predictive  technique. 

A  more  adaptable  method  of  simulating  the  viscous-inviscid  inter¬ 
actions  that  occur  in  typical  nozzle  flowfields  involves  solving  the 
time  dependent,  compressible  Navier-Stokes  equations  uniformly  over  the 
entire  nozzle  flowfield.  This  approach  has  a  direct  advantage  over  the 
previously  discussed  iteratively  patched  methods  where  an  accurate  vis¬ 
cous-inviscid  matching  procedure  is  required  in  order  to  obtain  reasonable 
results.  In  the  direct  approach,  the  predominantly  inviscid  and  viscous 
flow  regions  are  computed  simultaneously  with  no  matching  required. 

Holst  (5)  used  this  approach  to  solve  for  supersonic  flow  over  axisym- 
metric  boattail  nozzles  with  plume  simulators.  Although  a  plume 
simulator  does  not  model  the  entrainment  and  blockage  of  a  jet  plume, 
its  flowfield  does  contain  phenomena  characteristic  to  coflowing 
nozzles  such  as  turbulent  boundary  layers,  recompression  shock  waves, 
and  separated  recirculating  regions  of  flow.  Holst's  solutions  were 
obtained  using  MacCormack's  explicit  finite  difference  algorithm,  a 
stretched  mesh  aligned  with  the  solid  body  through  an  analytic  trans¬ 
formation  and  a  two  layer  eddy  viscosity  model  to  account  for  the 
Reynold's  stresses  that  included  a  relaxation  formula  to  model  the 
separated  flow  region.  Pressure  distributions,  skin  friction  co¬ 
efficients  and  areas  of  separated  flow  were  in  good  agreement  with 
experimental  data,  part icular i ly  in  the  cases  where  a  fine  mesh  was 
utilized.  Mikhail  (6)  recently  computed  solutions  for  viscous 
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supersonic  flow  around  an  axisymmetric  boattail  nozzle  with  a  jet 
exhaust  flow.  MacCormack's  explicit  method  was  again  used  as  the 
numerical  algorithm,  together  with  a  surface  oriented  mesh  system 
obtained  through  a  numerical  mapping  procedure.  Reynold’s  stresses 
were  also  accounted  for  through  the  application  of  algebraic  eddy 
viscosity  models.  Reasonable  agreement  with  experimental  surface 
pressure  data  on  the  boattail  was  obtained. 

Navier-Stokes  solutions  are  especially  useful  for  predicting  the 
off-design  nozzle  characteristics  where  the  flowfield  is  in  either  a 
significantly  overexpanded  or  underexpanded  state.  At  these  conditions 
the  flow  structure  is  usually  more  complex  with  viscous  regions  becoming 
more  prevalent  than  at  on-design  conditions.  One  feature  typical  of  these 
off-design  conditions  is  the  establishment  of  a  triple-point  in  the  jet 
flow,  and  the  appearance  of  a  strong  normal  shock  wave  referred  to  as  a 
Mach  disc  in  axisymmetric  flow  or  a  Riemann  wave  in  two-dimensional  flow 
(Figure  1).  This  strong  shock  formation  occurs  when  the  deflection  angle 
of  the  jet  flow  is  large  enough  so  that  the  resulting  shock  wave  is  too 
strong  for  a  regular  reflection  at  the  centerline  to  exist.  Near  the 
centerline  the  Mach  disc  must  be  normal,  since  this  is  the  only  way  a 
shock  can  occur  without  any  change  in  flow  direction.  As  shown  in  Figure 
2,  both  the  Mach  disc  and  the  reflected  shock  are  curved  near  the  triple 
point  (7).  A  slip  line  emanates  from  the  triple  point,  and  the  flow 
downstream  of  the  Mach  disc  and  reflected  shock  is  rotational  in  nature 
due  to  the  curvature  of  the  shocks.  As  discussed  by  Henderson  and  Lozzi 
(8),  this  region  downstream  of  the  shocks  may  be  either  totally  subsonic 
or  contain  both  supersonic  and  subsonic  regions.  If  the  incident  Mach 
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igure  1.  Shock  Structure  for  a  Typical  Overexpanded 
Coflowing  Axisyrmnetric  Nozzle. 
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REFLECTED 


ematic  of  a  Typical  Mach  Reflection 


number  is  greater  than  2.40,  region  3  will  be  supersonic,  while  region 
4  remains  subsonic.  Since  the  jet  Mach  numbers  relevant  to  this  investi¬ 
gation  are  greater  than  2.40,  the  latter  case  in  which  a  subsonic  core 
only  exists  behind  the  Mach  disc  will  be  examined. 

The  transition  from  regular  to  Mach  reflections  can  also  be  examined 
using  a  hodograph  diagram  shown  in  Figure  3.  For  deflection  angles  less 
than  the  transition  angle  (6  ) ,  the  flow  can  be  brought  to  the  required 

zero  deflection  at  state  3  through  a  weak  regular  wave  reflection.  Tor 
deflection  angles  greater  than  the  transition  angle,  the  flow  in  region 
3  cannot  achieve  a  zero  deflection  state,  and  lies  on  the  strong  shock 
portion  of  the  initial  shock  polar.  The  flow  near  the  centerline  passes 
through  a  strong  normal  shock  to  condition  4'  with  no  deflection  occurr¬ 
ing.  The  flow  state  at  the  curved  Mach  disc  then  exists  along  the  strong 
shock  portion  of  the  incident  shock  polar  from  a  zero  deflection  state 
4'  near  the  centerline  to  a  deflected  state  4  near  the  triple  point  with 
a  pressure  and  flow  angle  equal  to  that  in  region  3,  but  with  different 
velocity  and  entropy  values  that  generate  the  slip  line. 

The  mixed  supersonic-subsonic  flow  region  surrounding  the  Mach 
disc  greatly  complicates  the  anlysis  of  nozzle  flowfields  in  which  this 
shock  structure  is  present.  Flowfields  containing  this  phenomena  have 
not  been  adequately  simulated  in  the  past  using  viscous  techniques. 

A  variety  of  techniques  for  locating  the  triple  point  and  the  result¬ 
ing  normal  shock  have  been  presented  that  utilize  an  iterative  combina¬ 
tion  of  the  method  of  characteristics  and  schemes  involving  approximate 
analyses  such  as  pressure  requirements  downstream  of  the  strong  shock 
(9,10)  or  one-dimensional  flow  calculations  downstream  through  a  throat 
region  in  the  flow  (11,12).  Although  each  of  these  methods  give  reason- 
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Figure  3.  Hodograph  Diagram  Showing  Shock  Polar  Intersections 
!  for  Transition  from  Regular  to  Mach  Reflection. 

) 
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able  results  for  determining  the  triple  point  location  and  size  of  the 
resulting  Mach  disc,  each  is  only  valid  for  a  limited  range  of  nozzle 
pressure  ratios  and  jet  Mach  numbers.  In  addition,  none  of  these 
techniques  give  a  solution  for  the  subsonic  core  region  downstream 
of  the  normal  shock.  At  least  two  time  dependent  inviscid  techniques 
have  been  used  to  overcome  the  deficiencies  of  the  semi-empirical  methods 
previously  mentioned.  Jofre  (13)  performed  a  finite  difference  technique 
for  an  underexpanded  jet  with  a  Mach  disc  solution.  A  method  of  charac¬ 
teristics  solution  was  used  in  the  plume  expansion  region  near  the 
nozzle  exit.  This  gave  an  upstream  flow  profile  used  in  the  time  depend¬ 
ent  solution  further  downstream.  Sinha,  Zakkay,  and  Erdos  (14)  analyzed 
a  two-dimensional  underexpanded  jet  containing  a  strong  normal  shock 
using  a  finite  difference  technique  over  the  entire  flowfield  of  interest. 
Both  of  these  investigations  used  versions  of  Lax-Wendroff  numerical 
algorithms  and  simple  square  grids.  These  solutions  were  much  more 
adaptable  than  the  previous  semi-empirical  techniques  since  the  flow 
tends  to  adjust  to  its  local  environment  so  that  the  proper  shock  stucture 
is  automatically  obtained  as  the  solution  develops.  However,  these 
solutions  represent  only  a  first  approximation  of  the  correct  viscous 
solution  (15),  particularly  in  the  region  downstream  of  the  normal  shock 
as  shown  in  Figure  4.  Flow  properties  in  this  region  were  found  to  be 
heavily  dependent  on  the  level  of  damping  used  in  the  solution  algorithm. 
For  example,  Jofre  (13)  found  that  an  unrealistic  region  of  reverse  flow 
was  generated  immediately  behind  the  normal  shock  unless  heavy'  damping 
was  applied  in  the  solution  procedure. 

All  of  these  inviscid  solutions  involved  jets  exhausting  into  a 
quiescent  atmosphere.  Solutions  were  not  attempted  for  the  more  com- 
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plex  case  of  a  coflowing  nozzle  where  the  external  flow  stream  interacts 
with  the  jet.  A  full  Navier-Stokes  solution  which  accounts  for  the  viscous 
effects  present  in  the  flow  field  at  these  off-design  conditions  is 
necessary  in  order  to  adequately  simulate  both  the  strong-shock  structure 
with  its  resulting  imbedded  subsonic  flow  region  as  well  as  the  inter¬ 
action  of  the  jet  plume  with  the  external  flow  field. 

1 . 2  RESEARCH  OBJECTIVES 

The  primary  objective  of  this  research  is  the  development  of  a 
numerical  Navier-Stokes  method  capable  of  accurately  predicting  super¬ 
sonic  coflowing  nozzle  flowfields  which  contain  both  highly  viscous 
regions  and  complex  shock  structure  typified  by  the  Mach  disc  shock 
formation.  Overexpanded  axisymmetric  nozzles  will  primarily  be  simu¬ 
lated,  since  they  meet  the  previous  criteria  while  possessing  fairly 
compact  flow  domains  which  contain  the  flow  phenomena  of  interest.  The 
experimental  data  of  Bromm  and  O'Donnell  (16)  has  been  chosen  as  a  basis 
for  comparison  in  this  research  effort.  Data  in  this  reference  is  given 
for  an  axisymmetric  Mach  three  isentropic  nozzle  embedded  in  a  turbulent 
Mach  1.94  external  flowfield  as  shown  in  Figure  5.  Nozzle  pressure 
ratios  ranging  from  a  slightly  underexpanded  condition  to  a  highly 
overexpanded  condition  which  exhibits  the  Mach  disc  structure  were 
obtained  experimentally.  This  particular  nozzle  possesses  a  relatively 
thick  base  annulus  which  generates  a  strong  viscous-inviscid  interaction 
in  the  near  wake  region  of  the  nozzle.  These  interactions  affect  the 
development  of  the  primarily  inviscid  shock  structure,  and  can  only  be 
analyzed  properly  using  a  full  Navier-Stokes  methodology. 
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Although  Mikhail  achieved  full  Navier-Stokes  solutions  for  an 
axisvmmetric  coflowing  nozzle,  he  was  not  able  to  generate  an  accurate 
solution  for  the  condition  at  which  a  Mach  disc  shock  structure  was 
shown  to  exist  expe~imentally  (6).  Possible  causes  of  this  inability 
to  generate  the  strong  shock  structure  include  boundary  condition 
formulation,  mesh  spacing  and  turbulence  modeling.  These  three  areas 
will  be  concentrated  on  in  the  present  investigation  in  order  to  achieve 
the  desired  goal  of  an  accurate  predictive  technique  for  both  on-design 
and  off-design  nozzle  performance. 
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CHAPTER  II 


MATHEMATICAL  DESCRIPTION  OF  THE  FLOW  STRUCTURE 


2 . 1  GOVERN  INC  EQUATION’S 

The  governing  equations  for  flows  containing  the  shock  and  vis¬ 
cous  phenomena  of  interest  are  the  conservation  equations  for  mass, 
momentum,  and  energy  known  as  the  Navier-Stokes  equations.  The  gases 
involved  are  assumed  to  be  single  component,  have  constant  specific 
heats,  and  obey  the  perfect  gas  equation  of  state: 

P  =  pRT  (2-1) 

In  computational  fluid  dynamics,  the  Eulerian  method  is  usually 
applied  to  the  problem  of  interest.  This  method  involves  a  fixed  con¬ 
trol  volume  that  is  specified  relative  to  a  given  coordinate  system. 
Properties  of  the  fluid  are  then  specified  as  functions  of  both  time 
and  space.  The  conservation  equations  are  approached  using  this  method¬ 
ology. 


CONSERVATION  OF  MASS 


For  a  given  system  in  which  matter  is  neither  created  or  destroyed, 
the  law  of  mass  conservation  can  be  written  as 


pu)  dV  =  0 


(2-2) 


where  V  is  an  arbitrary  volume  fixed  in  space. 


CONSERVATION  OF  MOMENTUM 

For  a  given  system,  the  law  of  momentum  conservation  states  that 
the  rate  of  change  of  momentum  is  equal  to  the  sum  of  the  external  forces 


acting  on  the  control  volume.  If  body  forces  are  neglected,  this  law 
can  be  written  as: 


f/Mf 


(ou)  u ]  dV  = 


(V  •  S)  dV 


The  variable  S  denotes  a  stress  tensor  involving  pressure  and  viscous 
forces  which  acts  on  the  fluid. 


CONSERVATION  OF  ENERGY 


The  law  of  conservation  of  energy  states  that  for  a  given  system 
which  does  not  contain  any  internal  heat  sources,  the  rate  of  change 
of  the  total  energy  of  the  system  is  equal  to  the  heat  added  into  the 
system  plus  the  work  done  on  the  system  by  viscous  and  pressure  forces 
This  can  be  stated  as 


ffh 


(oe)  u] 


V  =  JJJ ( V«u *S  -  V* 


q)  dV 


Since  these  conservation  equations  are  valid  for  any  arbitrary 
volume  V;  when  the  integrands  are  continuous,  these  equations  imply  that: 


c* r  /  o t  +  v  •  pu 

=  0 

-»■ 

(2-5) 

d(,;u)/i'<t  +  V  • 

(cu) 

u  - 

v  •  s  =  0 

(2-6) 

d(.e)/dt  +  V  * 

(re) 

* 

u  + 

(V*q  -  V*u*S)  =  0 

(2-7) 

It  should  be  noted  that  these  equations  are  written  in  conserva¬ 
tive  form  where,  for  the  two-dimensional  and  axisymmetric  flows  of 
interest,  the  applicable  dependent  variables  are  p,  pu,  pv,  and  pe.  As 
showTi  by  Roache  (17),  this  conservative  form  allows  the  finite  difference 
equations  to  preserve  the  Gauss  divergence  property  of  the  continuum 
equations.  This  form  allows  a  balance  between  the  flux  quantities  and 
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accumulation  rates  for  a  small  control  volume.  Roache  also  states  that 
the  Rankine-Hugoniot  shock  relations  were  derived  using  the  conservative 
form.  Thus,  shock  jump  conditions  are  automatically  satisfied  since  the 
conservative  variables  are  continuous  across  the  shock  and  need  no  special 
treatment  because  of  d iscont init ies .  This  approach  is  known  as  shock 
capturing  or  shock  smearing.  The  conservation  form  of  the  equations  then 
allows  the  finite  difference  formulation  to  satisfy  the  physical  laws  on 
a  macroscopic  scale,  not  merely  in  some  academic  limit  as  Ax,  Ay,  and  At 
approach  zero. 

Since  the  flowfields  of  interest  are  turbulent,  the  solution  of  the 
conservation  equations  must  take  into  account  the  effects  of  the  random 
fluctuations  of  the  dependent  variables  inherent  to  turbulent  flows.  In 
accounting  for  these  effects,  cartesian  tensor  notation  will  be  applied. 
The  usual  conventions  of  a  repeated  subscript  indicating  summation  over 
the  entire  range  of  indices  and  a  comma  representing  partial  differentia¬ 
tion  will  be  used  to  make  the  equations  compact.  Cartesian  tensors  are 
used  to  allow  working  directly  with  the  physical  components,  while  still 
being  applicable  to  the  2-D  and  axisymnetric  systems  of  interest.  The 
conservation  equations  (2-5)  through  (2-7)  can  then  he  written  as: 


If. 


panded  into  the  following  form: 


In  these  expansions  the  barred  variables  represent  time  averaging  over 
a  time  interval  that  is  long  compared  to  turbulent  eddy  fluctuations, 
yet  small  compared  to  macroscopic  flow  changes.  The  primed  variables 
then  represent  fluctuations  due  to  the  turbulent  nature  of  the  flow. 

As  discussed  by  Chapman  (18),  this  time  averaging  approach  is  valid 
since  the  frequencies  of  most  unsteady  flows  of  interest  are  a  factor 
of  10  to  100  below  the  mean  frequency  of  turbulent  eddies. 

If  the  dependent  variables  u,  v,  and  e  are  mass  averaged  as 
described  in  reference  (19),  and  p  and  P  are  mean  (time  averaged)  state 
variables,  then  the  conservation  equations  can  be  written  in  the  form  of 
mean  flow  equations  as: 


c ,t  +  (pUj)  ,  =  0 

(2-13) 

(‘n.),t  +  K'u.u.)  +  Pi..  - 

(*ij  "  "  0 

(2-14) 

(re),  +  [peu,  +  q.  +  pule' 

t  J  J  J 

- s,  <>1:l  -  5-ixjn.j  - » 

(2-15) 

where  a  higher  order  mean  energy  dissipation  term  in  ul^  has  been  neglected 
in  the  energy  equation  (2-15). 

The  term  [-pujuj]  is  known  as  the  Reynolds  stress.  It  represents 
a  momentum  transfer  caused  by  turbulent  fluctuations  present  in  the  flow- 
field.  This  Reynolds  stress  term  can  be  written  as  an  apparent  stress 


caused  by  the  turbulent  nature  of  the  flow: 
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ij  turb 


-  . 1  u  .  u  . 
i  J 


(2-16) 


Since  air  is  essentially  an  isotropic  fluid,  the  mean  stress  term  can 
be  expanded  into  its  normal  and  shear  stress  components  as: 

(2-17) 

The  turbulent  stress  term  can  then  be  written  in  analogous  form  as: 


t  .  .  =  >u,  .  6 .  .  +  p  (u .  .  +  u .  . ) 

ij  k,k  ij  i,j  j,i 


1..  L  ,  =  A  u,  ,  6..  +  t  (u.  .  +  u.  .) 

lj'turb  t  k,k  ij  l,j  j,i 


(2-18) 


where  >.  and  t  are  the  turbulent  viscosity  coefficients  of  the  flow. 

The  coefficient  r.  is  known  as  the  eddy  viscosity,  and  is  analogous  to 
the  molecular  viscosity  coefficient  p .  However,  t  is  more  a  property  of 
the  dynamics  of  the  flow,  whereas  t  is  only  a  property  of  the  fluid. 
Combining  the  mean  and  turbulent  stress  terms,  an  overall  stress  term 
can  be  written  as: 


T  .  .  =  ( A+  >  )  u,  ,  5 .  .  +  ( U+£ )  (u .  .  +  U  .  .  ) 

IJ  !  total  t'  k,k  IJ  1,J  J  ,i 


(2-19) 


In  the  energy  equation  an  additional  unsteady  term  appears.  This 
term  is  by  nature  an  apparent  heat  flux  caused  by  the  fluctuations  in¬ 
herent  to  turbulent  flow  and  can  be  written  as: 

q.  L  u  =  cule’  (2-20) 

j  |  turb  j 

If  the  heat  flux  term  q  is  defined  by  the  Fourier  heat  equation 


q j  =  “kT. j  =  -  (Cpp/Pr)  T, ^ 


(2-21) 


then  by  the  former  analogy  q^ 


turb 


can  be  written  as: 


q . 


.  I  ,  ,  =  -  (C  f./Pr  )  T,  . 

J  turb  P  t  ’j 


(2-22) 


where  again  f  is  the  eddy  viscosity  coefficient,  and  Pr  is  the  turbulent 
Prandtl  number  of  the  flow.  Combining  these  two  heat  fluxes,  a  total 
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heat  flux  can  be  written  as: 


1j  total  *  -CP  (l'/Pr  +  T'j 


(2-23) 


The  mean  conservation  equations  can  then  be  written  in  the  follow¬ 
ing  form,  where  the  overbars  on  the  terms  are  dropped  for  convenience, 
and  where  the  values  of  the  shear  stresses  and  heat  fluxes  are  the  total 
values : 

P,t  +  (iu.)  =  0  (2-24) 


(;u.),t  +  [(cu.)u.  +  P6..  -  t..],.  =  0 
(oe),t  +  1  (re)  uj  +  '  ui  1 ij 5  -  j  =  0 


(2-25) 


ue),t  +  l  (re)  uj  +  -  ui  t  ]  =  0  (2-26) 

Since  the  flowfields  of  interest  are  either  two-dimensional  or 
axisvmmetr ic  in  nature,  the  mean  conservation  equations  can  be  written 
in  the  following  compact  vector  form: 

tir  l  "i  ( ^2 0  r  \  n 


HI  IF  1 _  3  (rJ  °G)  H 

HT  rjD  3r  J°  '  fjo 


(2-27) 


where  jc  =  0  or  1  for  either  the  two  dimensional  or  axisvmmetric  cases, 
respectively,  and 


r  1 

i 

| 

:  u  :  _ 

cu"  -  0 

;  f  = 

XX 

.  v  i 

CUV  -  T  : 

1 

,  xr 

r-e  , 

cue  +  q  -  uo  -  vt  i 

x  xx  xr  J 

(2-28) 


„  1 1:  uv  -  t 

G  "j  9  Xr 


iVe  +  q  -  ut  -  vo 

r  xr  rr, 


where 


=xx  =  -p  +  (>+>t>  div  V>  +  2^  +  t)  g 


t  =  -P  +  (-‘+*^)  div  V  +  2  (  u  +  t) 


(2-29) 

(2-30) 


(A  +  At)  =  -2/3  (u  +  £.)  (2-37) 

The  governing  equations  for  the  problem  of  interest  now  consist 
of  the  four  conservation  equations  in  matrix  form  (eq.  2-27)  with  four 
unknown  dependent  variables  p,  pu,  pv  and  pe.  The  perfect  gas  law  is 
used  to  define  the  pressure  in  terms  of  these  conservative  variables, 
and  a  model  of  the  dependence  of  the  eddy  viscosity  on  the  mean  flow  must 
be  introduced  to  overcome  the  "turbulent  closure"  problem  inherent  in  the 
turbulent  mean  conservation  equations. 

For  numerical  computation  in  a  transformed  (p,r,)  cartesian  plane,  the 
matrix  form  of  the  conservation  equations  (2-27)  can  be  written  as: 


,  '  (rJ  °G) 

r 


+ 


[  n 


X 


n 


r 


•■(rj"G)1 

Dn 


(2-38) 


where  r  and  ri  are  now  the  independent  variables,  and  the  transformation 
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derivatives  f.  ,  £  ,  n  ,  and  n  are  obtained  numerically  from  a  mapping 
procedure.  Equation  (2-28)  is  actually  in  weak  conservation  form  due 
to  the  varying  coefficients  in  front  of  the  derivatives,  and  also  due 
to  the  source  term  in  the  axisymmetric  case. 

2 . 2  BOUNDARY  AMD  INITIAL  COX'D  IT  ION’S 

Boundary  and  initial  conditions  must  be  given  in  order  to  solve 
the  conservation  equations  which  govern  the  flowfield.  These  conditions 
must  be  carefully  specified,  since  many  flow  features  such  as  shock  waves, 
boundary  layers,  and  recirculation  areas  arise  from  boundary  conditions. 
For  the  solution  of  a  symmetric  two-dimensional  or  axisymmetric  supersonic 
jet  embedded  in  a  supersonic  external  flowfield,  the  domain  of  interest 
can  be  defined  as  shown  in  Figure  6.  Only  one-half  of  the  total  nozzle 
flowfield  needs  to  be  considered  due  to  the  axis  of  symmetry  on  the  jet 
centerline.  The  remainder  of  this  chapter  will  detail  the  specific  bound¬ 
ary  conditions  that  are  pertinent  to  this  problem. 

THJE  _rPS_TR!wVr_  BOFNDARY 

Inflow  conditions  on  the  upstrean  jet  boundary  (AB)  and  the  up¬ 
stream  external  boundary  in  the  freestream  (CD)  are  completely  specified. 
Velocity,  pressure,  and  temperature  profiles  determined  from  auxiliary- 
computations  or  known  experimentally  on  this  boundary  fix  p,  pu,  fv, 
and  i-e  for  the  duration  of  the  problem  solution  at  the  boundary. 

THF_ FPPER  BOCNPARY 

The  upper  boundary  (DF.)  must  accurately  represent  a  free  flight 
condition  where  mass  flow  is  allowed  across  this  boundary  embedded  in 
the  supersonic  external  flowfield.  Weak  shock  waves  and  brand t 1 -Mover 
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gure  6.  Physical  Domain  for  the  Computational  Solu 


expansion  waves  must  also  exit  this  boundary  without  being  reflected 
artificially  back  into  the  domain  of  interest.  One  conditio’  that 
allovjs  this  is  the  assumption  of  a  simple  wave  solution: 

vr  =  0  (2-39) 

dAl'DE 

where  >  is  the  straight  left  running  characteristic  line  passing  through 
each  point  on  the  upper  boundary.  This  characteristic  line  is  determined 
only  by  the  value  of  the  Mach  number  and  flow  angle  of  the  supersonic 
flow  present  near  the  upper  boundary.  This  condition  assumes  that  the 
flow  along  this  boundary  is  inviscid  and  homentropic. 


THE  DOWNSTREAM  BOUNDARY 

The  downstream  boundary  (EF)  is  unique  in  that  no  rigorous  assumptions 
can  be  made  about  either  the  variables  or  their  gradients  unless  the  bound¬ 
ary  is  placed  a  great  distance  downstream.  In  this  case  a  no  change 
condit ion 


jf 


=  0 


(2-40) 


jx  |  EF 

could  be  assumed,  where  f  denotes  the  primitive  variables  p,  u,  v,  and 
T.  For  the  case  where  the  downstream  boundary  is  placed  where  gradients 
do  exist,  an  extrapolation  method  based  on  this  fact  can  be  reasonably 
applied.  One  such  method  is  to  assume  that  a  flow  gradient  accurate  to 
second  order  can  exist.  This  can  be  implemented  as: 


3  •  f  j 

C.x3— ^)|  =0  (2-41) 

jx  'EF 

In  other  words,  gradients  can  occur  which  are  parabolic  with  respect 
to  x.  This  condition  is  reasonable  if  the  gradients  at  this  boundary 
are  not  severe  as  in  the  case  where  a  strong  shock  wave  exits  the  boundary. 
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Tin:  r.r.NTr KLINE 


1 1 10  centerline  boundary  (AF)  is  a  line  of  symmetry  with  no  mass 
or  energy  flux  across  it.  Therefore,  the  following  boundary  conditions 
can  he  applied 


v  Iaf  =  ° 


(2-42) 


(2-43) 


(2-44) 


Since  the  v  component  of  velocity  is  zero  on  the  centerline,  this 
boundary  is  also  a  streamline  in  the  jet  flow.  For  steady,  adiabatic  flow 
with  negligible  volume  forces,  the  total  enthalpy  along  any  streamline 
is  a  constant.  Therefore,  along  the  centerline. 


Ho  |  =  (C  T  +  — )  =  constant 


(2-45) 


Since  the  condition  at  the  jet  exit  is  specified,  the  centerline  boundary 
can  be  properly  posed  using  this  approach. 


NOZZLE  WALLS 

The  nozzle  walls  (BG,  GH,  and  CH)  are  considered  to  be  no-slip, 


impermeable  boundaries.  This  assumption  gives  the  conditions  that: 


(2-46a) 


V I wa 1 1  =  0 


(2-46b) 


Since  the  stainless  steel  nozzle  consists  of  thin-walled  material 
with  a  thermal  conductivity  much  greater  than  that  of  the  surrounding 


fluid,  the  nozzle  walls  are  assumed  to  be  at  a  constant  temperature: 
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T  i  =  constant  (2-47) 

|wall 

This  wall  temperature  is  determined  by  applying  a  heat  flux  bal¬ 
ance  across  the  jet  and  freestream  boundary  layers  as  outlined  in 
Appendix  A. 

The  pressure  on  each  nozzle  wall  is  unknown,  but  can  be  approxi¬ 
mated  by  applying  the  degenerate  form  of  the  appropriate  normal 
momentum  equation  at  each  nozzle  surface  to  obtain  the  following: 


c?n 

In  this  expression  n  is  the  direction  normal 
s  is  the  direction  parallel  to  the  surface. 

INITIAL  CONDITIONS 


cU 


sn 


wa  1 1 


wall 


(2-48) 

to  the  wall  surface. 


and 


Since  the  governing  equations  contain  time  dependent  terms,  initial 
conditions  must  be  specified  before  the  solution  process  can  begin.  The 
specification  of  these  initial  conditions  is  somewhat  arbritrary  in  nature, 
althougn  steep  gradients  must  be  avoided  to  prevent  numerical  divergence 
during  the  solution  process.  Since  the  flow  is  predominantly  supersonic 
in  nature,  the  incoming  flow  profiles  will  have  a  dominant  effect  on  the 
solution  in  the  whole  computational  domain.  The  incoming  flow  conditions 
are  imposed  on  the  complete  domain  as  discussed  in  the  section  on  initial 
condition  implementation  of  chapter  5. 
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CHAPTER  Ill 


NUMERICAL  PROCEDURE 

The  nu.ner ical  procedure  consists  of  solving  the  governing  equa¬ 
tions  with  applicable  boundary  and  initial  conditions  through  the  use 
of  appropriate  finite  difference  techniques  on  a  high-speed  computer. 

This  procedure  can  be  broken  down  into  several  elements  which  include 
the  finite  difference  coordinate  system,  the  solution  algorithm,  and 
the  convergence  criteria  used  in  the  computational  solution.  Each  of 
these  areas  will  be  discussed  in  this  chapter. 

3  .  1  tTMlRDIflATE  SYSTEM 

Dqma in  of  Computation 

The  physical  domain  of  computation  consists  of  a  rectangular  area 
defined  by  orthogonal  coordinates  (x,r)  as  shown  in  Figure  7.  The  mesh 
consists  of  II  points  on  the  x  axis  and  JL  points  on  the  r  axis,  where 
II.  and  .11.  are  dependent  on  the  extent  of  the  physical  domain  required  for 
the  particular  jet  plume  case  of  interest. 

Mesh  Stretching 

In  order  to  obtain  an  accurate  numerical  solution  of  a  viscous  flow- 
field,  the  mesh  spacing  must  be  much  finer  in  areas  containing  relatively 
high  gradients  of  the  variable  properties  such  as  velocity,  density,  and 
temperature.  In  the  coflowing  nozzle  these  high  gradient  areas  include 
the  boundary  layers  on  the  nozzle  walls  and  the  shear  layer  in  the  wake  of 
the  nozzle  annulus.  This  stretching  is  accomplished  through  the  use  of  a 
patched  exponential  stretching  scheme  of  the  following  form: 


r(j)  =  I. 

in 


Cr,<  j  ) 

(e _ -H 

(eC-  1) 


f or  j 


(3-1) 


where  I.  and  r,  are  defined  as  shown  in  Figure  8.  The  constant  C  is 

in 

determined  bv  the  minimum  spacing  Ar  .  desired  for  the  mesh  next  to 

min 

the  wall  boundary.  Applying  the  desired  A r  .  to  equation  (3-1)  gives 


Ar 


C  = 


*7  ( 1  + 

on  e 


jniri 

L 


(e 


1,) 


(3-2) 


m 

The  value  of  C  is  then  obtained  through  the  use  of  an  iterative 
Aitken  extrapolation  technique. 

This  mesh  stretching  procedure  is  applied  in  the  radial  direction 
on  both  sides  of  the  nozzle  wall  where  boundary  layers  are  present.  It 
is  also  applied  in  tlu  axial  direction  at  the  end  of  the  nozzle  where 
the  jet  flow  begins  to  expand  or  contract  and  where  the  near  wake  due 
to  the  nozzle  annulus  begins  to  form. 


Adaptive  Mesh 

It  is  desireable  that  the  fine  merh  remain  in  the  areas  of  rela¬ 
tively  high  velocity  and  temperature  gradients  as  the  solution  progresses 
towards  convergence.  This  is  not  a  problem  in  the  case  of  boundary 
layers  that  are  adjacent  to  a  fixed  wall,  but  is  a  concern  in  the  free 
shear  layer  area  generated  by  the  nuzzle  wake  and  the  interaction  between 
the  jet  and  freestream  flows.  This  shear  region  on  the  jet  plume  boundary- 
will  deflect  to  a  degree  that  is  primarily  dependent  on  the  nozzle  pressure 
ratio.  The  fine  mesh  region  should  therefore  also  be  adapted  to  conform 
to  the  general  position  of  the  shear  region. 

Hirt  (21)  has  used  a  technique  in  the  solution  of  free  surface  flows 
that  allows  the  grid  to  adapt  as  the  solution  progresses.  The  following 
kinematic  equation  is  applied  in  the  region  where  the  shear  layer  is  present. 


Figure  8.  Fxponentially  Stretched  Mesh  Schemati 


3r 

3t 


_  ,  3r. 

C  (v  -  u  — ) 
A  dx 


(3-3) 


This  equation  insures  the  condition  that  as  the  solution  converges 
the  physical  slope  of  the  constant  n  finite  difference  cell  boundaries 
is  the  same  as  that  of  the  velocity  vectors  near  each  cell.  When 
applied  in  a  finite  difference  format,  the  grid  can  then  adapt  to  the 
placement  of  the  shear  layer  as  shown  in  Figure  9.  Details  of  this 
process  are  explained  in  Appendix  B. 


Coordinate  Transformation 

The  physical  domain  as  typically  shown  in  Figure  9  is  mapped  to  a 
unit  square  in  the  computational  plane  shown  in  Figure  10.  The  con¬ 
stant  n  lines  are  aligned  parallel  to  the  centerline  and  the  constant 
f,  lines  are  in  the  direction  normal  to  the  centerline.  The  numerical 
algorithm  operates  on  this  coordinate  system  using  the  transformed 
conservation  equations  (2-32).  Care  must  be  taken  in  generating  the 
physical  mesh  so  that  smoothness  of  the  transformation  coefficients 
(£  ,  £r,  n  ,  and  n  )  is  retained  in  order  to  reduce  numerical  errors 
caused  by  the  mesh  configuration. 

3.2  SOLUTION  ALGORITHM 


MacCormack's  Method 

The  weak  conservative  form  of  the  two-dimensional,  time-dependent 
Navier-Stokes  Equations  (eqn.  2-38)  is  solved  using  MacCormack’s  ex¬ 
plicit  finite  difference  method  (22).  This  algorithm  is  an  efficient 
Lax-Wendroff  type  differencing  scheme  of  second  order  accuracy  which 
utilizes  time-splitting  and  two  step  predictor-corrector  techniques. 
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MucCormack's  algorithm  was  chosen  for  application  to  the  nozzle  problem 
because  of  its  previous  success  in  computing  inviscid-viscous  interact¬ 
ing  flows,  its  stability  in  supersonic  flow,  and  its  computational 
efficiency  achieved  by  time-splitting  the  finite  difference  operators. 

The  computational  solution  is  advanced  in  time  by  applying  the 
numerical  operator  to  the  solution  of  the  flowfield  at  time  t.  This 
can  be  written  as: 

I’U,  t  +  At)  =  I. (At )  •  uu,  n,  t)  (3-4) 

where  L(At)  is  the  two-dimensional  numerical  operator  representing 
Macformack ' s  algorithm  acting  on  the  transformed  conservation  equa¬ 
tions.  Through  the  use  of  a  time-splitting  technique,  this  two- 
dimensional  operator  l.(At)  is  separated  into  two  one-dimensional  sweep 
operators  in  the  t.  and  n  directions.  The  operator  !, .  (At)  denotes 


the  solution  of  the  equation: 


+  •  4  +  -4-  r  =  0 


r 


in  the  ■'  direction  hv  time  increment  At  seconds.  Similarly,  the 
operator  1.  (At)  represents  tilt  solution  of 


,  1,1  :-(rJoG) 

—  +  .  -  •  +  — —  n - : - 

•t  X  "  Jo  r  dr. 


in  the  Ti  direction  by  a  time  increment  At  seconds.  The  dependent 
variable  vector  !'(■',  r,,t)  can  then  be  advanced  in  time  as 


n  t  +  ;.t)  =  [i 


.''  /2(.M/M) -L  ( .',t )  •  L. M  / 2 ( At  /M)  ]  •!'(£,  n,  t 


)  (3-7) 


.t  =  At.  if  At.  ■  At 
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or  as 
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where  V1?  .  is  the  known  value  at  time  t,  and  l'n+ is  the  intermediate 
1 .J  i,J 

predictor  value.  The  actual  computed  value  at  time  (t+At)  is  then 
calculated  by  applying  the  following  corrector  algorithm; 


r  .  =  1/2  {  r .  .  +  vr;  -  (— )  i(5  ).  .  (fvTt  -.  -  f.  f> 
i.j  i.j  i,j  L-  y-  i»j  i+i, x  i,j 


+  4-  (?.)..  (rj°  .  Cn+J  .  -  rj°.  g7+:-)  1  > 

rJ°  r  l.j  1+1,  j  1+1,  J  i,J  i ,  J 

i.j 


(3-10) 


In  this  £.  sweep  predictor-corrector  algorithm,  the  matrices  F  and 
G  are  functions  of  the  following  difference  quotients: 


pm  _  p1*1  _  p™ 

rm  rm  =  /  /rm  i+l,i  i.j  i , j+1 

*i,j  ’  ’i,j  Ui,j’  U.  ’  2An  ' 


(3-11) 


This  scheme  can  be  depicted  graphically  as  shown  in  Fig.  11. 

The  L  (At)  numerical  sweep  operator  is  formulated  in  a  similar 
manner.  The  intermediate  predictor  value  is  given  by  the  expression: 


un+:  =  rn  .  -  (£4  t(n  ).  .  (Fn  .  -  rn  .  .) 

i ,  j  i ,  j  An  xi, j  i , j  i , j -1 

+  -4-  (n  ).  .  (4°.  Gn  .  -  4°.  Gn  .  .)]  + 

Jo  r  1,J  1,J  i,j  1.J-1  1,J-1 


jo  At  H 


(3-12) 


The  corrector  value  at  time  (t+At)  is  then  given  by: 

i’n+1  =  1/2  l  r"  .  +  rf’:  -  44  [(r,  ).  .  (Fn+/x1  -  Fn+4 
i.J  i.J  i,j  An  x  i,j  i,j+l  i,j 


„n+  i 
J  o  At  H  .  . 


+  ~  (n  ).  .  (4°  G .  j  .  -  rJo.  Gn  :)  ]  + - -  --’-J  )  (3-13) 

rJo  r  i,j  i,J+l  i,J+l  i,J  i.J  fJo 

i  ,  j 

The  matrices  F  and  G  are  now  functions  of  the  following  difference 


quotients : 


Fm  cm  .  =  /(r".1 

i ,  j  i ,  j  i ,  j 


rm  ,  .  -  rm 


i’m  ...  -  um 


1+1  « j _ LtL-1  i  .  j  +  1  j  ,  j 


(3-14) 
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£  SWEEP 


PREDICTOR  CORRECTOR 


7?  SWEEP 

PREDICTOR  CORRECTOR 


pure  11.  Graphical  Representation  of  the  Numerical 
Sweep  Operators. 


This  scheme  is  also  depicted  graphically  in  Figure  11. 


Ttie  KacCormack  algorithm  satisfies  the  requirements  of  stability, 
consistency,  and  accuracy  if  the  following  conditions  are  met:  It  is 
table  if  the  time  step  limit  satisfies  the  CFL  stability  condition  given 
in  the  next  section.  It  is  consistent  if  the  sum  of  the  time  steps  for 
each  of  the  sweep  operators  is  equal,  and  the  algorithm  is  second  order 
accurate  if  the  sequence  of  sweep  operators  is  symmetric. 


and 


At  -  1/2  - 

r  '  1 

o  vPr  Pr  ; 


(3-18) 


Finally,  analysis  of  the  mixed  derivative  terms  found  in  the  Navier- 
Stokes  equations  gives  the  condition: 


At  =  At 
x  r  - 


2  x  A  r 


1 


(3-19) 


[-0  +  >t)  (v  +  O] 


For  the  finite  difference  equations  applied  to  the  complete  set  of 
Navier-Stokes  equations,  the  stability  criteria  can  then  be  estimated  a 

(3-20) 


At  •  minimum  (At  ,  At  ) 
-  .  .  x  r 

i  -J 


where 

At 

x  — 


i"!  s  <Fr*-fFt,+b  *  on 


(3-21) 


and 

At 


!v  i  +  c  +  |  t  ^  )  +  I-  [-(>+>t)(p  +  0]’J}  (3-22) 


For  the  present  non-Cartesian  jet  flow  cases  that  were  computed, 
the  maximum  time  step  was  calculated  as: 


At  =  minimum  (At . ,  At  ) 
'  n 

t  ,j 


where 
.At  _  = 


AS. 


u  .  1  +  c  +  •-  i  Tg  (jl~  +  p— )  +  I  -  (  3+ 1 1 )  ( P  +  £•)]  *  r 

t  n 


(3-23) 


(3-24) 


and 

At  = 


AS 


;  +  c  *  i  ■  <iv  +  pt  > +  Jr  Po+^Ki.  + 1 )  i  - 1 

ii  t  n 


(3-23) 
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0-26) 
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-  x .  .  )“  + 
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(3-27) 

n 

i  ,J 

i  ,J 

u  .  = 
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(u  .  .  ,  V  .  .  , 
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(3-28) 
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) 
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> 
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•H 

l.J-1 

(3-29) 

The  actual  time  step  used  in  the  numerical  procedure  was  less  than 
this  estimated  maximum.  Factors  varying  between  0.35  and  0.80  were  used 
during  the  computations.  Flow  solutions  involving  relatively  large  viscous 
wakes  containing  recirculation  regions  required  much  smaller  allowable  time- 
steps  with  CFL  factors  on  the  order  of  0.35-0.40. 

M'MF.P.  ICAI.  DAMPING 

Strong  shocks  imbedded  in  a  flowfield  being  solved  computationally 
can  often  cause  numerical  oscillations  (17)  which  may  lead  to  program 
failure  due  to  physically  unrealistic  values  of  computed  pressure, 
density  or  temperature.  These  oscillations  are  caused  by  numerical 
truncation  errors  and  can  be  reduced  by  refining  the  grid  in  the  areas 
of  shock  locations.  However,  this  can  be  impractical  when  the  oscilla¬ 
tions  are  of  a  transient  nature  caused  by  computational  start-up  or 
re-start  procedures,  or  where  the  shock  location  varies  for  different 
experimental  cases  and  mesh  refinement  for  each  individual  case  is 
undesireable .  In  this  situation,  a  fourth  order  pressure-gradient 
damping  concept  as  introdu'  .iCCormack  and  Baldwin  (23)  can  be 

applied  to  increase  the  stability  of  the  numerical  algorithm. 

This  damping  scheme  is  applied  in  both  the  £  and  n  directional 


37 


sweeps.  In  the  "  numerical  sweep,  damping  is  implemented  by  replacing 

F  .  bv  (F. .  .  +  F  )  and  G. .  .  by  (G. .  .  +  )  in  equation 

H>3  ii, J  D  n, J  ii,  3  D  /  M 

11  jJ  1  1  >  J 

(3-5).  The  predictor  and  corrector  steps  in  this  case  are  represented 
by  i i  =  i  and  ii  =  i+1,  respectively.  The  damping  terms  F^  and  are 
in  the  following  form: 

|P. .  -  2P.  .  +  P.  ,  .1 

n  =  a  (ju  |  +c  >  I  1+1 _ —La _ JdLtiL 

D.  .  .  (.  1  ii, j  1  ii,y  (P...  .  +  2P.  .  +  P.  ,  .) 

ii, J  i+l, J  i,J  1-1, J 


O'  !  .  -  U.  .) 

1+1, J  1,J 


(3-30) 


P.  .  .  . .-  2P. .  .  +  P. .  .  .  i 

li  >3+1 _  n,3  n.,1-1 


n  (|Vii,j'  +  Cii,j}  (P  +  2P..  .  +  P..  . 

11,3+1  n,3  n,3-l 


(U  -  U  )  .  ° 

'  ii , j+1  ii,j-l;  ii, j 


(3-31) 


In  the  n  sweep,  a  similar  procedure  is  implemented  by  replacing  F.  .. 

1  >  J  J 

by  (F.  . .  +  F  )  and  G.  . .  by  (G.  . .  +  G  )  where 
1,33  ,JJ  i  ,jj 


■X.  (  ju  .  .  .  +  C  .  .  .  )  4— 

r-  1,33  1  ,33  (P 


I  i+l, j j 


2P.  . .  +  P. 


. .  +  2P.  . .  +  P.  .  . . ) 
i+l,33  i,33  i-l,33 


i+l  ,33  Li-l,jj) 


(3-32) 


,P.  . . ,  -  2P.  .  +  P.  .  . 

1  1,3+1  i,3  i,3-l  I 


■;  =  a  (  v  +  c  )  J  "’4  '  ' - - x  »-j.  2LI 

D.  ..  n1  i  ,j j  '  i.jy  (P.  . .,  +  2P.  .  +  P.  .  . ) 
i,3J  i,3+l  i ,3  1,3-1 


(U .  . -  r.  .)  ‘  rJo . . 

1,3+1  i,3  i,33 


(3-33) 


In  this  case  the  predictor  and  corrector  steps  are  represented  by 
jj  =  j  and  jj  =  j+l,  respectively.  For  both  sweeps  a  and  u  are 
damping  constants  where  normally 


ir  =  a  =0.5 
5  h 


(3-34) 
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3 . 3  CONTERC.ENCE  CRITERIA 

Convergence,  as  applied  in  this  section,  refers  to  iteration  con¬ 
vergence  as  opposed  to  truncation  convergence,  which  involves  the 
convergence  of  the  solution  of  the  FDE  to  the  solution  of  the  PDE  as 
Ax,  Ar,  and  At  0.  Iteration  convergence  refers  to  the  arrival  at  a 
solution  to  the  discretized  h'avier-Stokes  equations  within  some  accept¬ 
able  tolerance  through  the  use  of  an  iterative  process.  As  stated 
by  Roache  (17),  there  are  no  definitive  criteria  for  iteration  con¬ 
vergence.  A  somewhat  subjective  judgement  of  convergence  must  be  made 
based  upon  an  examination  of  the  iterative  behavior  of  the  solution 
flow  variables.  Different  flow  variables,  as  well  as  variables  at 
different  locations,  converge  at  different  rates.  If  the  slowest 
converging  variable  in  the  flow-field  is  known,  it  should  be  the  most 
closely  examined  for  convergence. 

In  the  present  case  for  a  coflowing  supersonic  nozzle  with  a  re¬ 
latively  thick  base  annulus,  an  examination  of  the  flow  variables 
revealed  that  the  slowest  converging  variable  was  the  base  pressure 
of  the  nozzle  annulus.  The  location  of  this  base  pressure  is  within 
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the  subsonic  flow  area  involving  recirculation  in  the  wake  of  the  nozzle 


annulus  as  shown  in  figure  12.  The  flow  variables  in  this  subsonic 
region  converged  much  more  slowly  than  did  those  in  the  predominantly 
supersonic  jet  and  freestream  flows. 

Since  the  coflowing  nozzle  problem  primarily  involves  high 
Reynolds  number  flow,  the  advective  terms  in  the  conservation  equa¬ 
tions  dominate  the  viscous  diffusion  terms.  A  characteristic  time  for 
a  disturbance  to  cross  the  flowfield  may  then  be  characterized  by: 


'ch  u 


(3- 3b) 


ch 


where  L  is  the  length  of  the  flowfield  in  the  direction  parallel  to 
the  characteristic  velocity  u^.  For  the  jet  problem  u  ^  was  repre¬ 
sented  by  u ^ .  Since  in  general  the  magnitude  of  was  less  than  ujet> 
this  gave  a  more  conservative  estimate  of  the  characteristic  time. 

The  convergence  criteria  was  then  established  by  the  following 
procedure.  The  numerical  solution  was  either  initially  started,  or 
restarted  from  a  previous  case.  As  the  solution  converged,  the  base 
pressure  was  monitored  until  its  magnitude  varied  less  than  +  1%  for 
one  characteristic  time  period.  At  the  end  of  this  characteristic  time 
period  the  solution  was  stopped  as  convergence  was  achieved.  Visual 
comparison  of  Mach  Number  and  density  profiles  over  the  flowfield  con¬ 
firmed  the  convergence  of  the  solution  using  this  procedure.  A  sample 
base  pressure  convergence  plot  is  shown  in  figure  13. 
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M  MH1  H  Oh  mHATIo\S  X  10 


Figure  13.  Typical  Base  Pressure  Convergence,  P./P  =  0.527 


CHAPTER  IV 


BOUNDARY  AND  INITIAL  CONDITION  IMPLEMENTATION 

As  in  chapter  two,  boundary  and  initial  conditions  must  be 
defined  in  order  to  solve  the  conservation  equations  which  govern 
the  flowfield.  Values  of  the  dependent  variables  for  points  on 
the  boundaries  of  the  computational  domain  must  be  specified  in 
order  for  the  numerical  operators  to  solve  the  flowfield  correctly. 
This  section  presents  the  numerical  formulation  of  the  boundary  and 
initial  conditions  used  for  the  solution  of  the  coflowing  nozzle. 
The  conditions  were  presented  in  a  mathematical  context  in  chapter 
two . 

4 . 1  THE  UPSTREAM  BOUNDARY 

The  flow  properties  on  the  upstream  boundary  (AB  and  CD  of 
Figure  6)  are  held  fixed  for  the  duration  of  the  computational 
solution.  The  values  of  these  properties  were  derived  in  the 
following  manner.  In  the  external  flowfield,  a  parabolized  Navier- 
Stokes  solution  (24)  was  computed  for  the  ogive  body  used  in  the 
experimental  coflowing  nozzle  tests  as  shown  in  Figure  14.  This 
solution  determined  that  the  pressure  gradient  at  the  inflow  bound¬ 
ary  in  the  external  flow  stream  in  negligible,  and  that  the  static 
pressure  at  the  inflow  boundary  is  99%  of  that  in  the  undisturbed 
flow  in  the  wind  tunnel.  The  static  pressure  along  the  ogive  body 
surface  shown  in  Figure  15  was  then  used  as  an  input  to  a  two- 
dimensional  turbulent  boundary  layer  code  (25,  26)  along  with  the 
other  specified  freestream  conditions  (M^  T^,  Re^,  T  )  to  produce 
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ssurc  Variation  along  the  Ogive 
1  i  zed  Navi er-Stokes  Solver  (23). 


input  velocity  component  and  temperature  profiles  at  the  upstream 
boundary.  The  conservative  variables  on  this  boundary  (CD)  are  then 
calculated  using  these  profiles  and  the  static  pressure  along  the 
boundary.  Flow  variables  on  the  upstream  boundary  in  the  jet  flow 
(AB)  are  determined  in  a  similar  manner.  The  same  boundary  layer 
code  is  applied  using  the  jet  exit  conditions  and  the  length  from 
the  nozzle  throat  to  the  nozzle  exit  plane  as  the  boundary  layer 
starting  length.  Again,  profiles  for  the  velocity  components  and 
temperature  are  obtained  along  the  boundary.  Values  for  the  conser¬ 
vative  variables  are  then  computed  using  these  profiles  and  the  value 
of  the  pressure  at  the  jet  exit.  Since  the  value  of  the  vertical 
velocity  component  is  zero  on  the  centerline  boundary,  a  polynomial 
fit  is  used  to  set  the  vertical  velocity  profile  from  the  edge  of  the 
boundary  layer  to  the  centerline. 

4  .  2  TI!K_  1'PPER  BOUNDARY 

T!ie  upper  boundary,  labeled  DE  on  Figure  6,  utilizes  the  simple 
wave  procedure  outlined  by  Roache  (17,  pp  282-283).  This  procedure 
assumes  that  properties  are  constant  along  a  straight,  left-running 
characteristic  line  passing  through  each  point  on  the  upper  boundary. 
Tiie  position  of  this  line  running  through  a  boundary  point  (i,JL)  is 
determined  by  the  angle  (p j  +  C),  where 

=  arcsin  (1/M.  .)  (4-1) 

M  t 

is  the  local  Mach  angle  for  supersonic  flow,  and 

=  arc tan  (v/u)  (4-2) 

is  the  local  fiow  direction.  The  properties  on  this  characteristic 
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line  are  determined  by  linear  interpolation  involving  the  properties 


at  points  (  i-1  ,  ,1b),  (i-1,  Jb-1),  and  (i,  .11.-1)  as  shown  in  Figure 
16.  Points  (i-1,  .11.-1)  and  (i,  .11.-1)  are  points  interior  to  the 
computational  domain,  and  point  (i-1,  II.)  is  either  a  point  on  the 
inflow  boundary  (i  =  1)  or  a  previously  defined  point  on  the  upper 
boundary  resulting  from  the  left  to  right  sweep  along  the  boundary 
using  this  technique. 

As  shown  in  Figure  16,  one  of  two  possible  interpolation  schemes 
is  applied  depending  on  the  local  values  of  the  quantities  ( f.  +  ') 


and  (Ar/Ax).  For  the  case  of  tai.  (u.,  +  ).  ,  ,,  , 

M  i-l,  JI.-l 


Ar/Ax,  the 


position  of  the  characteristic  line  lies  between  the  points  (i-1,  JL-1) 
and  (i,  JI.-l).  The  properties  at  the  point  p,  and  thus  at  point 
(i,  Jl.)  ,  can  then  be  determined  by: 

f  i  ,.11.  =  fp  =  f  i-1, JL-1  +  (Ix)(fi,JL-l  ”  f i-1, JL-13  (4-3) 

The  value  of  '•  ,  and  thus  the  position  of  the  point  p,  can  be  deter¬ 
mined  by  the  following  procedure.  If  the  quantity  w  is  defined  by: 

w  =  tan  [90°  -  ( ux,  +  )]  ( A— 4 ) 

then  by  geometry: 

w^  =  (Ax  -  v)/Ar  ( 4  —  5) 

If  the  interpolation  procedure  of  equation  (4-3)  is  applied,  then 

w  =  w  +  (— )  (w.  -  w.  )  (4-6) 

p  i-l, JL-1  ox  i,.JL-l  i-l, JL-1 

Equating  (4-5)  to  (4-6)  and  solving  for  >  gives : 


(Ax/ Ar)  -  w.  .  IT  , 

_ _ i-l  , JL- 1 _ 

(w.  -  w  .  )  / A.x  +  ( 1  / A. r  ) 

l  , JL-1  i-l , JL-1 


<  (.'.r /Lx)  the 


For  the  alternate  case  of  tan  ( ;.  + 

in  i-l,JL-l 

position  of  the  characteristic  line  lies  between  the  points  ( i — 1 
Jl.-ll  and  (i-1,  II.)  .  The  properties  at  p  are  then  determined  by 
equat  ion : 

fi,JL  fp  fi-l, .11,-1.  +  /“r)  (f i-1  ,JL  "  fi-l,.JL-l) 

In  this  case,  the  quantity  w  is  defined  as 

w  =  tan  (;.  +  ) 

m 

and  by  geometry 

w  =  c.r  -  ;  )/:.x 
p 

Tile  interpolation  scheme  for  w  now  gives  the  expression: 

WP  =  Wi-1„IL-1  +  ('/:-r)  (wi-l,JL  -  Wi-1,JL-1) 

Again,  equating  (4-10)  to  (4-11)  and  solving  for  >  gives: 

(Ar/ix)  -  vi_i)jI_1 


(w.  .  -  w.  .  „  .  )/Ar  +  (1/Ax) 

t-1, JL  l- 1 , JL- 1 


the 


(4-8) 


(4-9) 


(4-10) 


(4-11) 


(4-12) 


Applying  the  computed  values  of  the  interpolation  length  to  the 
respective  interpolation  equation  (either  (4-3)  or  (4-8))  gives  the 
proper  value  of  the  desired  flow  variable  at  the  boundary  point  (i, 
JL)  . 


4  .  3  THE  DOV.’hSTREAM  BOUNDARY 

The  downstream  boundary  (EF)  is  placed  in  a  region  where  gradients 
in  the  flow  variables  are  expected  to  exist.  A  quadratic  extrapula- 
tion  can  be  used  on  this  boundary  that  lets  3f /3x  and  \<~  {  /  be 
nonzero,  thus  satisfying  this  gradient  condition.  Assuming  a  constant 
grid  spacing  near  this  boundary,  Taylor  series  expansions  can  be 
performed  in  the  following  manner  from  a  point  (II,.  j)  on  the  houndarv 
to  the-  following  three  points  interior  to  the  boundary: 


*  IL-l.j  =  f  11.,  j  +  'X  *x  !  1 1. ,  j  + 


o 

■*x^ 


'3  4  I 

TI  1  +  0(Z-X  4t,  •>  r4-13> 

IL’J  4x3!il’j 


f  T1.-2 ,  j  ■  r„.,j  *  =•*  %Ui +  2“»2  r'liiuj +  <’<4*5  -*%.,>  «*«> 

-  /v  d  X 


o  2 

,  'f  j  9 Ax"  f 

f  ,,  .  =  f , .  .  +  i.'.x  .  + 

1L-3,j  11.,  j  ;  I L ,  j 


■JX 


2  !  I L ,  j 


3  '4  I 

,  +  C(fx  -4)  L,  .)  (4-15) 


3x 


3  !  IL.  j  ’ 


If  the  assumption  is  made  that  the  last  term  in  each  equation  can 

be  neg 1 ec t ed ,  i . e . 

3  4 

Ax  =  0  (4-16) 


then  equations  (4-13)  through  (4-15)  can  be  solved  simultaneously  to 
give  the  following  expression  for  the  boundary  point  (IL,  j)  in  terms 
of  the  interior  points: 


fTI  .  =  3  fTt  .  .  -  3  f_,  ,  .  +  f  TT  ,  . 
1L,j  IL-l,j  IL-2.J  1L-3, j 


(4-17) 


This  condition  works  well  as  long  as  large  pressure  gradients 

do  not  exist  at  this  boundary,  as  in  the  case  when  a  normal  shock  wave 

3  3 

exits  the  boundary.  If  this  does  occur,  the  term  involving  3  f / 3x 
is  no  longer  neplegible.  Equation  (4-17)  then  can  become  numerically 
unst  ab 1 e . 


Numerical  divergence  did  occur  when  the  previous  extrapolation 
condition  was  applied  to  regions  of  subsonic  flow  present  at  this  out¬ 
flow  boundary.  Therefore,  the  following  first  order,  zero  gradient 
condition  was  applied  at  points  (IL,  j)  when  the  Mach  number  at  points 
(II.-l,  j)  vas  found  to  be  subsonic: 

f  t  i  •  =  fTT  ,  •  (4-18 

T  . .  ,  j  ■  1  1.  1  ,  j 

4.4  THE  CFTJFJU.IKE 

The  centerline  boundary  (AF  in  Figure  6)  is  a  line  of  svmmetrv 
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with  no  mass  or  energy  flux  across  it.  The  vertical  velocity  component 


condition  (eqn  2-36)  is  applied  by  setting: 

v.  =  0  for  1  •  i  '  IL  (4-19) 

-1,1  - 

The  symmetry  conditions  for  the  density  and  u  component  of  velocity 
(eqns  2-37  and  2-38)  which  are  valid  only  at  the  centerline  are  applied 
in  the  following  manner.  Taylor  series  expansions  are  performed  for 
.  and  u  from  a  centerline  boundary  point  (i,l)  to  points  (i,2)  and 
(i,3)  which  are  distances  Ar  and  KAr,  respectively,  above  the  center- 
lino  boundary.  The  series  expansions  to  these  points  give  the  following 
equations  for  density: 

i  0  3 

>i,3-  +  t!«.i  :%,»>  <4-20b) 

1  or  -r 

If  the  centerline  symmetry  condition  (2-37)  is  applied  and  the 
higher  order  term  in  each  equation  is  neglected,  these  two  equations 


can  be  solved  simultaneous!.'  to  obtain: 


i  ,1 


-  <■-'  1,2  -  ~  1,3 


(4-21a) 


(K“  -  1) 

A  similar  expression  can  be  obtained  for  the  u  component  of  velocity: 


u  .  -  u  .  0 

i,2  i,3 


(4-21b) 


ar  -  i) 

The  previous  extrapolation  boundary  conditions  for  the  density  and 
the  horizontal  velocity  component  were  applied  only  to  regions  of  super¬ 
sonic  flow  at  the  centerline.  t'ndesireable  pressure  wigcles  occurred  if 
these  conditions  were  applied  to  regions  of  subsonic  flow.  The  following 
first  order  zero  gradient  condition  was  applied  at  the  centerline  points 


(i,l)  when  the  Mach  number  at  points  (i,2)  was  found  to  be  subsonic 


*  i,l 


i  >  2 


(4-22a) 


u 


i,- 


( A  —  2  2  b ) 


Since  the  v  component  of  velocity  is  zero  on  the  centerline, 
this  boundary  can  be  considered  as  a  streamline.  As  discussed  in 
chapter  2,  the  stagnation  enthalpy  is  then  constant  on  this  boundary. 
Since  the  total  enthalpy  at  the  inflow  boundary  is  known,  then 


o  i,l 


H  .  .  for  1  <  i  <  IL 
o  1 , 1  — 


(4-23) 


Since  the  value  of  u ^  ^  has  been  previously  determined  by  equation 
(4-22),  the  definition  of  the  stagnation  enthalpy  can  be  expanded  to 
determine  the  value  of  temperature  at  each  boundary  point  (i,l): 


T.  =  T.  + 

i,l  1,1 


[  (u  )"  -  (u  )"] /2C 

1,1  i,l  p 


(4-24) 


The  values  of  the  primitive  variables  ,  ,  u,  v,  and  T  have  now 
been  determined  for  each  centerline  boundary  point,  so  that  the  re¬ 
quired  values  of  the  conservative  variables  can  be  computed  along 
this  boundary. 


4 . 5  THE  NOZZLE  WALLS 

The  nozzle  walls  (BG,  GH,  and  CH  in  Figure  6)  are  treated  as 
no-slip,  impermeable  boundaries.  The  no-slip  condition  is  imposed 
on  the  three  wall  faces  by  imposing  the*  following  conditions  (see 
F igure  17): 

Inner  wa 1 1 

u.  =  v.  =  0  for  1  •  i  <  IW  (4-25a) 

i,JUI  l , TW1 
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Outer  wall: 


u.  =  v.  =0  for  1  <  i  <  IW 

1 ,  JWO  l ,  JWO  — 


(4-25b) 


Base  (vertical)  wall: 


u  .  =  vTI,  .  =  0  for  JWI  <  j  <  JWO 
1W,J  IW,j 


(4-25c) 


As  discussed  in  chapter  2,  a  constant  wall  temperature  is  im¬ 
posed  on  the  nozzle  walls.  This  condition  is  applied  simply  as: 

h.JKI  '  Ti,JViO  '  TW  fpr  1  '  1  i™  <4'26a) 

and 

Tt1i  .  =  T. ,  for  JWI  <  j  <  JWO  (4-26b) 

IW,j  W 


A  first  order  pressure  gradient  condition  derived  from  equa¬ 
tion  (2-42)  is  applied  on  each  wall.  This  is  imposed  as: 


i ,  JWI 

=  Pi,JWI-l  for  1 

<  i  <  IW 

(4-27a) 

i  ,  JWO 

=  Pi , JWO+1  f0r  1 

<  i  <  IW 

(4—2  7b ) 

II 

•r-> 

PIW+l,j  for  JWI 

•  j  *■  JWO 

(4-27c) 

Since  the  points  (IW,  JWI)  and  (IW,  JWO)  are  positioned  at 
sharp  corner  poi-ts,  the  simple  pressure  condition  applied  in 
equation  (4-27)  is  i  applicable.  An  averaging  scheme  was  there¬ 
fore  used  to  allow  the  are  to  adjust  at  the  corners.  This 

averaging  is  applied  as: 


JWI 
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T!ie  primitive  variables  u,  v,  P,  and  T  have  now  been  defined 
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on  the  wo  1  1  boundaries.  1  he  required  values  of  the  conservative 
variables  can  then  be  computed  on  this  boundary. 


4.6  INITIAL  CONI) I TIONS 

As  discussed  in  chapter  2,  initial  values  of  the  conservative 
variables  must  be  imposed  over  the  computational  domain.  Since  the 
incoming  flow  variables  are  fixed  in  time,  they  are  initially  imposed 


over  the  complete  computational  domain  as  follows: 

V°.  .  =  l’°  .  ,  j  •  .JV.’I  and  JVO  j 
i.J  1,J  - 


14-24) 


The  value  of  the  u  component  of  velocity  in  the  wake  region  (JV.’I  ■  j 
JVO)  would  be  zero  from  the  input  profile.  Therefore,  the  u  com¬ 
ponent  in  the  wake  is  set  to  grow  exponentially  using  the  following 


equat ion : 


u.  =  k,,  (1  -  e  IV,j  i,j  B  ) 


where 
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for  IV  ■  i  •'  II.  and  JVI  •_  j  <_  JVO.  I'se  of  this  scheme  allows  the 
velocity  in  the  far  wake  to  be  close  to  that  of  the  two  streams,  thus 
acce 1 erat ing  convergence . 

Once  an  initial  case  of  the  coflowing  nozzle  had  numerically 
converged  to  a  valid  solution,  each  succeeding  case  was’  in  it ial ized 
by  imposing  a  new  jet  input  profile  on  the  inflow  boundary  of  the 
preceding  converged  solution.  This  technique  allowed  the  new  solu¬ 
tion  to  converge  at  a  much  greater  rate,  since  the  subsonic  recircu¬ 


lation  zone  in  the  near  wake  was  already  in  existence  f roc  the  previous 


sol ut  ion  . 


CHAPTER  V 


T l IRBI'LEXCE  MODELING 

The  experimental  tests  used  as  a  basis  for  the  computational 
solutions  were  conducted  at  a  Reynolds  number  of  2.2  x  10^  ,  based 
on  the  ogive  body  length  and  the  external  flow  conditions.  The 
external  flow  in  the  region  of  the  nozzle  is  thus  expected  to  be  of 
a  fully  turbulent  nature.  Reynolds  numbers  in  the  interior  jet  flow 
covered  a  range  from  1  x  10^  to  1.7  x  10^,  based  on  nozzle  exit  con¬ 
ditions  and  the  nozzle  throat  to  exit  plane  length.  Considering  the 
effects  of  compressiblity  and  the  existence  of  a  favorable  pressure 
gradient  in  the  divergent  portion  of  the  nozzle,  a  transition  Reynolds 
number  of  5  x  10^  was  used  to  determine  the  condition  at  which  the 
jet  flow  possessed  turbulent  characteristics  (27). 

The  turbulent  nature  characteristic  of  these  flows  can  be  accounted 
for  in  the  computational  solution  by  a  variety  of  eddy  viscosity  models 
ranging  from  locally  dependent  algebraic  models  to  the  more  complex 
higher-order  closure  models  such  as  the  turbulent  kinetic  energy  methods. 
Although  the  higher  order  methods  can  account  for  the  time  history  of 
the  turbulence  in  a  flow,  they  require  that  accurate  initial  profiles 
of  the  turbulent  shear  stress  be  known  or  reliably  calculated  (28) 

Tf  this  initiai  profile  condition  cannot  be  satisfied,  then  this  type  of 
prediction  method  cannot  be  effectively  utilized.  Since  this  proved  to 
be  the  case  for  the  jet  problem  under  consideration,  locally  dependent 
eddy  viscosity  models  were  carefully  applied  over  the  computational 
domain  . 
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As  shown  in  Figure  18,  the  computational  domain  contains  three 
distinct  regions  in  which  various  eddy  viscosity  models  are  applied. 
These  regions  consist  of  an  area  containing  boundary  layers,  a  far  wake 
region  downstream  of  the  nozzle  exit,  and  a  near  wake  region  close  to 
the  nozzle  exit  plane. 

5  • 1  BOUNDARY  LAYER  MODEL 

In  the  first  region,  the  dominating  flow  features  are  the  boundary 
layers  along  the  nozzle  walls.  Since  the  experimental  boundary  layer 
thicknesses  are  at  least  an  order  of  magnitude  smaller  than  the  nozzle 
radius,  a  two-dimensional  turbulence  model  was  judged  to  be  sufficient 
for  the  axisymmetric  cases.  The  eddy  viscosity  model  applied  in  this 
region  is  a  two  layer  Cebeci-Smith  model  (29).  The  inner  layer  of 
this  model  accounts  primarily  for  the  laminar  sublayer  adjacent  to 
the  wall,  with  the  outer  layer  accounting  for  the  remainder  of  the 
boundary  layer  region. 

The  expression  for  the  inner  model  is  based  on  Prandtl's  mixing- 
length  theory,  which  can  be  written  as: 


£  . 
i 


du 


(5-1) 


where  u^  is  the  local  tangential  velocity  parallel  to  the  wall  surface, 
and  r^  is  the  normal  distance  measured  from  the  wall.  The  mixing  length 
in  this  model  is  adapted  from  Van  Priest's  sublayer  model,  and  is  ex¬ 
pressed  as: 

-  r  ►  p^t 
n  w 

C.  =  0.4  r  (1  -  e  26,‘  )  (5-2) 

i  n 
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Figure  18.  I ddy  Viscosity  Model  Domains. 


€ 


Figure  19.  Kddy  Viscosity  Distribution  Across  a 
Boundary  Layer. 


58 


This  inner  model  assumes  the  flat  plate  pressure  condition  (dp/dx  = 

0) ,  but  could  be  modified  to  account  for  a  pressure  gradient  in  the 

direction  parallel  to  the  wall  within  the  sublayer. 

The  outer  region  eddy  viscosity  model  consists  of  a  Clauser-type 

displacement  thickness  model  defined  by  the  equation: 

cQ  =  0.0168  du  5  >.  ( 5-  3 ) 

e  1 

where  u^  is  the  appropriate  tangential  velocity  at  the  boundarv- layer 

edge  and 

j.  6  u 

t"  =  /  (1  -  — )  dr  (5-4) 

o  u  n 

e 

is  the  incompressible  displacement  thickness.  This  model  also  includes 

Klebanoff's  intermit tency  factor  defined  by  the  following  equation 

=  [1  +  5.5.  (y/O6!'1  (5-5) 

The  inner  and  outer  regions  of  each  boundary  layer  are  defined  by 

the  requirement  that  the  eddy  viscosity  remain  continuous  across  the 

entire  layer.  This  is  accomplished  by  applying  the  inner  model  outward 

from  the  wall  until  c.  =  cc  at  a  value  r  .  The  outer  model  is  then  applied 

l  c 

from  r  outward  across  the  remainder  of  the  flowfie.  i  the  boundarv 
c 

layer  region.  Figure  19  shows  a  typical  eddy  '•  v-.sit_  rofile  across 
this  region. 

5 . 2  FAPl  WAKE  MODEL 

In  the  region  downstream  of  the  nozzle  exit,  the  initial  boundary 
layers  on  the  nozzle  walls  merge  to  form  a  shear  layer  containing  an 
imbedded  wake  region.  This  region  in  the  flovfield  can  be  further 
divided  into  two  separate  areas:  the  near  wake  region  close  to  the 
nozzle  exit  that  contains  flow  features  such  as  the  corner  expansions, 
a  "deadwater"  zone,  recompression  stiocks,  and  the  far  wake  region  further 
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downst  ream  where  the  flow  eventually  tends  to  a  similar  free  shear  layer 


The  cutoff  value  used  to  define  a  representative  edge  of  the  mixing 
layer  was: 

.  ,  1  =  1  x  10A  sec-1  (5-10) 

edge 

Tli is  value  gave  a  reasonable  value  of  6  as  shown  by  a  typical  vorticity 
profile  in  Figure  20.  The  absolute  value  of  vorticity  typically  dropped 
to  less  than  one  percent  of  |u  i  within  one  cridpoint  outside  of  the 
cutof f  point  . 

A  flat  plate  validation  case  was  computed  using  the  far  wake  model 
to  check  its  accuracy  in  a  known  turbulent  flowfield.  The  data  of 
Toyoda  and  Hiriyarna  (33)  for  a  flat  plate  in  turbulent  flow  at  a  Mach 
number  of  1.6  was  used  as  the  basis  for  a  computational  solution.  Tire 
velocity  defect  in  the  wake  obtained  both  experimentally  and  numerically 
is  shown  in  Figure  21.  As  shown  by  this  figure,  the  results  generated 
by  the  numerical  turbulence  model  compared  very  well  with  the  thin  flat 
plate  data.  Further  details  of  this  computation  are  listed  in  Appendix 
C. 

5  . F.A R_KAKF  MODEL 

The  accuracy  of  a  Prandtl  type  mixing  length  turbulence  model  is 
substantially  dependent  on  the  use  of  length  scales  that  are  truly  re¬ 
presentative  of  the  flow  in  a  given  region.  In  the  far  wake  region 
where  a  single  mixing  layer  exists,  the  previous  definition  of  the  length 

scale  involving  :  is  valid.  However,  for  the  case  of  a  coflowing  nozzle 

w 

with  a  thick  base  annulus,  several  length  scales  need  to  be  defined  in 
the  near  wake  region.  As  shown  in  Figure  22,  the  existence  of  the  sub¬ 
sonic,  recirculating  "dead  water"  region  adjacent  to  the  nozzle  base 
wall  complicates  the  flow  simulation.  In  the  near  wake,  the  length 
must  transition  from  the  appropriate  boundary  laver  thickness 
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at  the  nott.lt  exit  plane  to  the  single  mixing  layer  thickness  (f  )  that 
exists  in  the  fa-  wake  region. 

Ihis  transition  is  accomplished  using  the  following  procedure.  The 
exterior  edges  of  the  mixing  layers  in  the  near  wake  are  determined  using 
the  vortieitv  limits  previously  defined  for  the  far  wake.  The  interior 
edges  of  the  dual  shear  layers  are  then  defined  by  the  Mach  0.5  contour 
line  surrounding  the  "dead  water"  region  as  shown  in  Figure  22.  The 
vortieitv  turbulence  model  defined  by  equations  5-6  and  5-S  is  then  applied 
in  the  near  wake  region  from  the  end  of  the  boundary  layer  zone  to  the  start 
of  the  initial  mixing  zone  with 

■  =  d  in  the  free  stream  flow 

o 

=  i .  in  the  jet  stream  flow 
l 

c  =  0.5  (c  +  i.)  inside  the  Mach  0.5  contour  of  the  "dead  water 
o  i 

region 

An  initial  mixing  zone,  one  base  height  in  length,  is  used  to  smoothly 

adiust  the  thickness  of  the  dual  shear  layers  (c  and  1.1  that  exist  in 

'01 

the  expansion  and  recompression  zones  of  the  near  wake  region  to  that  of 
the  single  shear  laver  (i  )  in  existence  further  downstream.  The  follow- 
ing  exponential  equation  is  applied  in  this  region: 


-k. 


Uy.)  =  d 


-  (i. 


o  or  j  x . 
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where  k^  is  determined  by  the  expression 

k.  =  2 (~—  -  -  -  0.5) 

1  yB  yB 


(5-11) 


( 5-12) 


and  the  value  of  x  lies  in  the  following  range: 

+  0.5v 


-  0.5y  •  x 


The  point  x  is  centered  in  the  mixing  region.  It 
the  position  of  the  midpoint  must  he  specified. 


determined  that 


If  a  1 lewed  to  11  oa t  , 


marine  r  . 


Figure  23.  Computed  Nozzle  Base  Pressure  vs  the  Position 
of  the  vixing  Region  Midpoint  x^. 
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CHAPTER  VI 


NTMERI CAL  RESULTS 

This  chapter  will  discuss  the  numerical  results  of  the  computa¬ 
tional  solutions  using  the  algorithm,  boundary  conditions,  and 
turbulence  modeling  detailed  in  previous  chapters.  The  first  section 
of  this  chapter  will  duscuss  the  experimental  cases  taken  as  the  basis 
for  comparison  with  the  numerical  solutions  for  the  coflowing  nozzle. 
The  next  section  discusses  details  involved  in  the  actual  computa¬ 
tional  procedure.  The  last  section  covers  the  comparison  between  the 
experimental  and  computational  solutions,  including  some  analyses 
of  the  accuracy  of  the  simulations  and  discrepancies  between  the 
numerical  solutions  and  the  experimental  data. 


6 . 1  Experimental  Data  Base 

As  outlined  in  Chapter  1,  the  data  of  Bromra  and  O'Honnell  (16) 
is  used  as  the  experimental  basis  for  this  research  effort.  Super¬ 
sonic  fields  of  flow  generated  experimentally  contain  both  highly 
viscous  flow  regions  as  well  as  shock  structures  ranging  from  weak 
regularly  reflected  shock  waves  to  the  strong  Mach  disc  shock  forma¬ 
tion.  Five  different  experimental  nozzle  pressure  ratio  conditions 
are  used  as  the  basis  for  the  computational  solutions.  Nozzle  base 
pressure  measurements  and  schlieren  photographs  are  the  basis  for 
experimental  versus  computational  comparisons. 

Mode  1 

The  model,  a  stainless  steel  body  of  revolution,  consisted  of 
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a  16.25  inch  ogive  forebody  and  a  one  inch  diameter  cylindrical 
afterbody  (see  Figure  15).  The  total  length  of  the  model  was  7.5 
inches.  The  model  was  supported  by  a  107  thick  hollow  side  strut 
that  acted  as  a  conduit  for  the  air  flow  to  the  jet  as  shown  in 
Figure  5.  The  effects  of  this  strut  were  found  to  be  negligible 
on  the  flow  field  downstream.  This  model  was  fitted  with  a  nozzle 
which  gave  essentially  isentropic  flow  with  an  exit  Mach  number  of 
3. 00.  The  inner  diameter  of  this  nozzle  at  the  exit  plane  was  0.742 
inches,  and  the  length  from  the  nozzle  throat  to  the  exit  plane  was 
1.20  inches.  Four  base  pressure  orifices  were  used  to  obtain  the  base 
pressure  measurements  as  shown  on  Figure  5. 

F.xternal  Flow 

The  external  flow  conditions  were  generated  in  the  NASA  Langley 
Q-inch  supersonic  wind  tunnel.  The  free  stream  Mach  number  was  set 
at  1.94,  and  the  free  stream  Reynolds  number  was  fixed  at  2.2  x  10^ 
based  on  the  body  length  of  the  model.  A  turbulent  boundary  layer 
on  the  model  at  these  conditions  was  insured  through  the  use  of  a 
transition  strip  near  the  nose.  The  tests  were  conducted  at  a 
tunnel  stagnation  pressure  of  one  atmosphere  (assumed  to  be  2116  psf). 
I'sing  these  given  conditions,  the  stagnation  temperature  of  the  free 
stream  was  calculated  to  be  equal  to  580. 5°R. 

Jet_  Flow 

The  flow  just  upstream  of  the  jet  exit  plane  was  given  to  be 
at  a  Mach  3,  zero  divergence  angle  condition.  The  total  pressure 
in  the  jet  flow  was  varied  to  obtain  the  desired  nozzle  static  pressure 
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ratio.  The  jet  static  pressure  at  the  nozzle  exit  plane  was  not 
measured  directly,  but  was  calculated  using  the  given  nozzle  area 
ratio  and  the  jet  total  pressure.  A  total  temperature  in  the  jet 
flow  was  not  given  experimentally,  but  was  assumed  to  be  equal  to 
the  freestream  stagnation  temperature  (5S0.5°R). 


6 . 2  Computational  Retails 

Solutions  were  computed  for  the  coflowing  nozzle  at  the  follow¬ 
ing  five  nozzle  static  pressure  ratios:  P^/P^  =  0.150,  0.251,  0.527, 
1.03  and  1.59.  These  solutions  were  all  performed  on  a  CDC  Cyber  175 
digital  computer  located  at  Wright-Patterson  AFB,  Ohio.  The  average 
rate  of  data  processing  was  0.0015  second  per  grid  point  per  iterative 
time  increment.  In  this  section  further  details  of  these  computations 
will  he  discussed. 


Cr_id  Parameters 

Important  parameters  of  the  computational  grid  including  the 
number  of  grid  points  and  axis  length  utilized  in  both  the  axial 
and  radial  directions  are  listed  in  Table  1  for  each  nozzle  pressure 
ratio  condition.  The  value  of  the  axial  field  length  includes  a 


length  increment  of  0.4  upstream  of  the  nozzle  exit  plane,  with  the 
exit  plane  at  a  value  of  x/r^  t  =  0.0.  A  compact  45  x  45  point  gric 
was  used  in  the  two  most  highly  overexpanded  cases  (P  /P^  =  0.150, 
0.251).  Twelve  additional  grid  points  were  then  added  in  the  axial 


direction  downstream  of  the  original  grid  to  form  a  57  x  45  grid  for 


the  next  case  where  P  /Vn  =  0.527  (Figure  25).  This  methodology  gave 
a  consistent  cell  length  in  the  axial  direction  for  each  of  these 
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three  cases  in  which  both  a  regular  shock  reflection  (P./P  =  0.527) 

J 

and  Mach  discs  (P  /P^  =  0.251,  0.150)  should  occur.  The  computational 
domain  was  then  stretched  using  the  57  x  45  point  grid  to  the  degree 
necessary  to  cover  the  phenomena  of  interest  for  the  two  underexpanded 


cases  . 

Minimum  grid  spacing  in  both  the  axial  and  radial  directions 
occured  adjacent  to  the  nozzle  walls,  and  was  set  at: 

Ax/r.  =  Ar/r.  =  0.030 
let  jet 

The  patched  exponential  stretching  outlined  in  chapter  3  was  then 
applied  to  form  each  grid  initially  for  the  various  cases.  The  adap¬ 
tive  grid  procedure  was  applied  in  the  radial  direction  during  the 
solution  procedure  as  discussed  in  Chapter  3  and  Appendix  B  to  obtain 
the  final  grid  geometry  for  each  case. 

Table  1.  Computational  Grid  Parameters 

P./P  IL  JL  XLT/r.  YLT/r. 

— J _ _ _ _ _ _ _ _ _ _ jet 


0.150 

45 

45 

3.4 

3.0 

0.251 

45 

45 

3.4 

3.0 

0.52  7 

57 

45 

4.4 

3.0 

1.03 

57 

45 

6.4 

3.0 

1  .  59 

57 

45 

8.4 

4.0 

Coarse  Grid  Effects  on  Boundary  Layer  Resolution 
The  numerical  solution  of  the  Navier-Stokes  equations  can  in¬ 
volve  significant  truncation  errors  in  regions  containing  high  velocitv 
gradients  such  as  within  turbulent  boundary  layers  when  relatively 
course  computational  grids  are  employed.  Errors  in  computed  velocity 
gradients  invol  ed  in  shear  force  terms  at  the  wall  can  result  in  erron¬ 
eous  values  of  the  pressure  gradient  along  the  wall.  The  extent  of 


these 


a  tur- 


s  can  be  realized  by  assuror,  ing  the  existence  of 
hulent  boundary  layer  possessing  a  velocity  profile  in  the  following 
form  (lb): 


+  + 
u  =  v 


u+  =  2.50  1 n ( y+ )  +5.10  11 
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(b. la) 


(6.1b) 


Since  the  dominant  term  in  the  shear  stress  is  the  gradient  , 

’V 

a  comparison  of  the  value  obtained  using  this  profile  can  he  made 
with  that  obtained  using  the  numerical  algorithm.  MacCortnack's  algorithm 
computes  this  gradient  using  a  first  order  finite  difference  in  the 
direction  of  the  sweep  for  the  predictor  step.  At  the  wall  this  grad¬ 
ient  is  computed  as: 


t  =  u  .  .  .  /  Ay 

y  '  w  i  ,  j w+1  w 


(b.2) 


where  jw+1  is  the  first  grid  point  above  the  wall,  and  Ay  is  the 


grid  spacin'.-,  adjacent  to  the  wall.  If  a  nominal  Reynolds  number  of 


1x10  based  on  body  length  is  applied  to  the  previous  turbulent  velocity 
profile,  the  ratio  of  actual  wall  velocity  gradient  given  by  the  pre¬ 
vious  profile  to  that  computed  using  equation  6-2  can  be  displayed  as 
a  function  of  the  ratio  of  u  velocity  component  at  the  first  point 
above  the  wall  to  the  freest  ream  velocity  as  shown  in  Figure  26.  For 
grid  spacing  substantially  greater  than  the  sublayer  (uj/u  =  0.4751, 
excessive  error  in  the  computed  velocity  gradient  occurs.  Since  in  the 
coflowing  jet  solutions  the  previously  stated  minimum  grid  spacing 

gives  values  of  (u./u  )  =  0.72  and  (u, /u  ) , .  =  0.90,  a  correction 

1  e  1  e  jet 

is  needed.  This  is  accomplished  by  applying  the  value  of  f  given  hv  a 
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( d u /d y ) | w  / Uu/ay)! 


a  hound a r v  layer  analysis  al  the  inflow  boundary  along  the  nozzle  wall 


when  j  =  jw,  instead  of  using  the  standard  differencing  procedure  in 
the  algorithm.  This  procedure  has  resulted  in  smooth  pressure  profile 
near  the  inflow  boundary  instead  of  slight  pressure  jumps  previously 
observed  in  course  mesh  cases.  This  concept  is  analogous  to  the  wall 
functions  described  bv  Launder  and  Spaulding  (37)  and  employed  bv 
Pier v  and  Forester  (38). 

External  Flow  Parameters 

Values  for  the  external  flow  variables  were  input  numerically 
bv  arriving  the  experimental  values  of  M  ,  P  ,  T  ,  and  the  boundarv 
layer  profiles  generated  for  use  on  the  inflow  boundary.  A  value  for 
t lie  turbulent  skin  friction  coefficient  of  C^.  =  0.00283  was  obtained 
from  the  boundary  layer  input  profiles  and  was  applied  on  the  ex¬ 
ternal  nozzle  wall  as  discussed  in  the  last  section.  A  wall  temper¬ 
ature  of  551 °R  was  calculated  using  the  analysis  in  Appendix  A. 

The  wall  temperature  was  given  a  constant  value  for  all  five  cases 
since  the  analysis  showed  that  this  temperature  should  not  vary  more 
than  one  degree  over  the  range  of  simulated  flow  conditions. 

let  !'  1  ow_  Paramo t vrs 

The  numerical  jet  field  of  flow  was  d'-tern.ind  bv  applying  the 
experimental  values  of  the  jet  Mach  number  and  total  temperature 
given  in  tin:  previous  section  along  with  the  values  of  the  jot  total 
pressure  and  skin  friction  coeflirient  listed  in  Table  2.  Since  tin 
nozzle  wall  tenperat  ure  was  considered  to  he  constant,  the  value 
for  the  internal  nozzle  wall  was  also  set  at  351  R.  As  shown  in  t  he 
t  a)  le,  the  state  of  the  jet  boundarv  layer  in  all  of  the  cases  t  >:- 
r  cjit  the  1  ov.-st  nozzle  pressure  ratio  condition  tp./p  -  n.lhni 


was  considered  to  be  turbulent  at  the  nozzle  exit  plane.  Since  the 


case  where  P./P  =  0.251  possessed  a  Reynolds  number  below  the  pre¬ 
viously  set  transition  point,  a  s  lotion  using  a  laminar  jet  boundary 


layer  was  also  obtained.  At  this  pressure  condition,  both  the  laminar 


and  turbulent  jet  boundary  layers  gave  nearly  identical  values  for  the 
nozzle  base  pressure  and  normal  shock  position.  However,  the  laminar 
case  exhibited  a  mild  shear  laver  oscillation  for  the  duration  of 


t  .e-  solution  procedure.  For  this  reason  the  results  of  the  turbulent 
case  are  presented.  It  is  interesting  to  note  that  no  mixing  layer 
oscillations  were  observed  in  the  other  laminar  jet  condition  (P./P  = 

0. 1  5:  i  . 
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significant  portion  of  the  flow  along  each  boundary  was  subsonic  in 
nature.  This  condition  only  occurred  for  the  highly  overexpanded 
case  (p7Pk  =  0.150),  where  a  substantial  embedded  subsonic  region 
exists  downstream  of  the  Mach  disc  shock  structure. 

The  quadratic  extrapolation  given  by  equation  4-17  produced  very 
reasonable  results  for  those  cases  where  the  outflow  was  either  entirely 
supersonic,  or  only  subsonic  at  a  very  few  points  in  the  nozzle  wake. 

For  the  case  where  a  substantial  area  of  subsonic  flow  existed  behind 
a  Mach  disc  structure  and  extended  to  the  donwstream  boundary,  num¬ 
erical  divergence  occurred  when  equation  4-17  was  applied.  A  second 
order  zero  gradient  condition  was  then  applied  in  regions  of  subsonic 
flow  along  this  boundary.  Application  of  this  condition  did  not 
produce  numerical  divergence,  but  did  give  unrealistic  pressure  jumps 
at  this  boundary.  A  first  order  zero  gradient  condition  given  in 
equation  4-18  was  then  successfully  applied  to  subsonic  regions  on 
this  boundary  with  reasonable  results. 

An  almost  identical  situation  occurred  along  the  centerline 
boundary  for  subsonic  regions  containing  fairly  strong  radial  flow 
gradients  close  to  the  centerline.  Both  the  extrapolation  condi¬ 
tion  given  by  equation  4-21  and  a  second  order  zero  gradient 
condition  produced  unrealistic  radial  oscillations  in  the  numerical 
solution  (called  ’’wiggles")  within  these  regions  of  subsonic  flow. 

The  application  of  a  first  order  zero  gradient  condition  given  by 
equation  4-22  helped  reduce  these  oscillations  to  achieve  a  reason¬ 
able  solution. 

Two  solutions  were  also  computed  where  first  the  downstream 
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boundary,  and  subsequently,  the  upper  boundary  were  repositioned 
greater  distances  from  the  nozzle  as  discussed  in  Appendix  E.  No 
changes  were  detected  in  either  the  shock  structure  or  the  nozzle 
base  pressure  coefficient.  These  test  cases  validate  the  effectiveness 
of  these  boundary  conditions  at  the  positions  utilized  in  the  actual 
nozzle  solutions. 

Convergence 

As  stated  in  chapter  3,  the  numerical  solutions  were  either 
initially  started  using  only  the  boundary  layer  profiles  across  the 
computational  domain,  or  restarted  from  a  previous  solution  by  apply¬ 
ing  new  input  profiles  at  the  jet  inflow  boundary.  Solution  times 
based  upon  the  convergence  criteria  discussed  in  chapter  three  varied 
significantly  for  the  two  methods  of  initial  startup.  Solution  times 
on  the  Cyber  175  using  only  the  boundary  layer  profiles  to  form  the 
initial  conditions  were  approximately  3.0  hours  for  the  45  x  45  point 
mesh  and  3.8  hours  for  the  57  x  45  point  mesh.  Solution  times  for 
cases  restarted  from  previous  solutions  were  approximately  1.7  hours 
for  the  45  x  45  mesh  and  2.1  hours  for  the  57  x  45  mesh.  The  large 
difference  in  these  solution  times  is  mainly  attributed  to  the  length 
of  time  required  for  the  subsonic  recirculation  region  in  the  near 
wake  to  form  and  achieve  a  steady  state  condition. 

6 . 3  Comparison  with  Experimental  Data 

Comparisons  between  the  numerical  solutions  and  the  experimental 
data  can  be  made  both  qualitatively  and  quantitatively.  Figures  27 
through  31  give  a  good  visual  comparison  between  the  numerical  solutions 
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depicted  as  Mach  number  contours  and  the  experimental  schlieren 
photographs.  In  these  figures  the  computed  solutions  above  the 
centerline  were  reflected  to  give  a  total  nozzle  flowfield  to  compare 
with  the  schlieren  photographs.  All  features  typical  of  afterbody 
types  of  flows  such  as  the  shock  structure  internal  to  the  jet  core 
flow,  external  recompression  shocks,  and  shaar  layer  development  are 
readily  discernible  and  in  very  good  agreement  with  the  experimental 
data . 

As  shown  by  the  previous  figures,  the  pressure  condition  at  which 

the  jet  flow  shock  structure  transitions  from  a  regular  reflection  on 

the  centerline  to  the  Mach  disc  formation  lies  between  the  two  cases 

where  P . /P  =  0.527  and  P./P  =  0.251.  Although  the  shock  structure 
J  oo  J  03 

near  the  centerline  appears  similar  in  the  computational  solutions 
for  these  two  cases,  an  enlargement  of  this  region  as  shown  in  Figure 
32  reveals  several  differences.  The  shock  strength  (related  to  the 
Mach  number  jump  across  the  shock)  is  much  greater  in  both  of  the  strong 
Mach  disc  cases  (P^/P^  =  0.150  and  0.251)  than  in  the  regularly  re¬ 
flected  case  (Pj/P^  =  0.527).  The  sonic  lines  in  this  region  are 
displayed  as  dashed  lines  in  Figure  32  in  order  to  easily  identify 
regions  of  subsonic  flow.  Both  Mach  disc  cases  contain  areas  of  sub¬ 
sonic  flow  downstream  of  the  shock  along  the  centerline,  whereas  the 
minimum  Mach  number  behind  the  regularly  reflected  shock  is  approxi¬ 
mately  equal  to  1.65. 

A  check  was  made  on  the  solution  to  the  highly  overexpanded  case 
(Pj /P^  =  0.150)  to  determine  if  the  numerical  solution  correctly 
simulated  the  flow  conditions  across  the  strong  normal  shock  in  the 
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COMPITKD  MACH  NUMBFR  CONTOURS  SCHLTFRFN  PHOTOGRAPH 


Figure  27.  Axisymmetric  Coflowing  Nozzle  Solution,  P./P  =  0.150 


COMPUTED  MACH  NUMBER  CONTOURS  SCHLIEREN  PHOTOGRAPH 


Axisymmetric  Coflowing 


region  near  the  centerline.  Since  the  v  velocity  components  are 
very  small  near  the  centerline,  a  one-dimensional  analysis  based 
on  the  Rankine-Hugoniot  relations  across  normal  shocks  can  be  applied. 
As  shown  in  Table  3,  the  computational  solution  was  within  2%  of 
the  exact  one-dimensional  analysis. 


Table  3.  Comparison  Between  a  1-D  Analysis  and  the  Computational 

Solution  Across  the  Mach  Disc  for  P./P  =  0.150. 

J  “ 


Exact  (1-D)  3.00  .475  10.33  3.857  2.679 

Computational  3.00  .442  10.52  3.906  2.695 

%  Error  1.1  1.8  1.3  0.6 

Several  other  phenomena  associated  with  afterbody  flows  are 
evident  in  Figure  33,  which  displays  computed  velocity  profiles 
at  given  axial  stations  for  the  large  Mach  disc  case.  The  separated 
"deadwater  zone"  of  recirculating  flow  is  readily  apparent  in  the 
near  wake  region,  as  is  the  development  of  the  near  wake  to  a  far 
wake  velocity  profile.  The  existence  of  the  strong  Mach  disc  near 
the  centerline  is  very  evident,  and  the  flow  in  the  subsonic  core 


region  behind  the  shock  accelerates  in  the  correct  manner  to  a  slight¬ 
ly  supersonic  condition  at  the  outflow  boundary. 

A  closer  look  at  the  near  wake  region  is  shown  in  Figure  34  for 
two  of  the  computational  cases.  This  figure  illustrates  the  change 
in  the  shape  of  the  "deadwater  region"  from  a  predominantly  symmetric 


nature  at  P./P  =  1.03  to  one  with  an  asymmetric  nature  at  P./P  - 

J  00  J  00 

0.150.  In  this  figure  the  dashed  lines  denote  the  dividing  stream¬ 


line  and  the  streamlines  through  the  stagnation  point  in  the  near 


wake  flow  for  each  case.  As  shown  in  Figure  35,  the  dividing  streamline 
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Figure  83.  Computed  Velocity  Profiles 


reamlinc  in  the  Compu 
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moves  toward  the  inner  wall  of  the  nozzle  as  the  total  pressure  in 
the  jet  is  decreased.  Although  the  stagnation  point  in  the  near  wake 
region  moves  radially  as  the  jet  stagnation  pressure  is  changed,  it 
remains  in  a  relatively  constant  position  axially  for  the  cases  com¬ 
puted  . 

Quantitative  comparisons  between  the  experimental  data  and  the 
computed  solutions  are  based  primarily  on  two  parameters:  the  axial 
distance  along  the  centerline  from  the  nozzle  exit  plane  to  the  point 
of  reflection  of  the  incident  shock  wave  at  the  line  of  symmetry, 
and  the  value  of  the  nozzle  base  pressure  coefficient.  This  reflec¬ 
tion  length,  along  with  the  type  of  shock  reflection  (either  strong 
or  weak) ,  is  a  good  indication  that  the  inviscid  flow  features  in 
the  jet  core  caused  by  viscous-inviscid  interaction  are  properly 
simulated.  These  computed  shock  reflection  lengths  are  obtained 
by  examining  the  axial  variation  in  Mach  number  along  the  centerline 
as  shown  in  Figure  38  for  a  typical  case.  As  shown  in  this  figure, 
the  shock  reflection  is  diffused  over  three  cell  lengths,  with  the 
computational  value  of  the  reflection  length  taken  as  being  at  the 
midpoint  of  these  three  cells.  Comparisons  between  the  experimental 
and  computational  values  of  these  reflection  lengths  are  shown  in 
Figure  37  and  Table  4.  Excellent  agreement  was  obtained,  with  the 
computational  results  being  within  2%  of  the  experimental  data. 

An  additional  quantitative  comparison  was  also  made  of  the 
computed  and  observed  Mach  disc  radii  for  the  two  cases  at  which 
the  Mach  disc  was  observed.  This  comparison  is  listed  in  Table  5. 

The  computed  Mach  radius  was  taken  as  the  radial  height  of  the 
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Figure  36.  Axial  Variation  in  Computed  Centerline  Mach  Number,  P./P^  =  0.150. 


sonic  line  immediately  behind  the  strong  shock  wave.  Very  good 


I 


agreement  was  obtained  for  the  highly  overexpanded  case  (P  /P  = 

3  00 

0.150)  in  which  a  fairly  large  Mach  disc  occurs.  Poorer  agreement 

was  obtained  for  the  case  that  was  much  nearer  to  transition  to  a 

regular  shock  reflection,  with  a  very  small  Mach  disc  radius  (r  /r. 

m  jet 

0.17).  The  first  cell  height  adjacent  to  the  centerline  in  the 
computational  solution  possessed  a  value  of  r/ret  ~  0.0 6,  so  that 
numerical  truncation  error  played  a  large  part  in  the  discrepancy 
in  Mach  stem  radius  at  this  particular  condition. 


Table  4.  Comparison  of  Shock  Reflection  Lengths 

P./P^  Experimental  Computational  %  Error 

3  x  /r.  (+0.05)  x  / r.  (+0.04) 

_ _ r  jet  — _ r  jet  —  _ 


0.150 

1.19 

1.17 

1 

O 

OJ 

0.251 

1.88 

1.91 

+0.5 

0.527 

2.98 

2.91 

-1.2 

1.030 

A. 37 

A. 26 

-1.9 

1.590 

5.6A 

5.53 

-2.0 

Table  5.  Comparison 

of  Mach  Disc  Radii 

P./P 

J  00 

Experimental 

r  /r.  (+0.05) 

m  jet  — 

Computational 

r  /r.  (+0.04) 

m  jet  - 

%  Error 

0.150  0.A5  0.AA  -2.2 

0.251  0.17  0.06  -2A . 
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Table  6.  Comparison  of  Base  Pressure  Coefficients 


P./P 

1  “• 

Exper imenta 1  * 
).'b(+0.003) 

Computat ional 

J.'B  (+0.002) 

%  Error 

0.150 

-0.240 

-0.310 

-25.0 

0.251 

-0.280 

-0.299 

-  6.8 

0.527 

-0.265 

-0.276 

-  3.9 

1.030 

-0.237 

-0.245 

-  2.9 

1.590 

-0.226 

-0.219 

+  2.5 

*interpolated  from  existing  values 

Comparisons  between  the  experimental  and  computational  values 
of  nozzle  base  pressure  coefficient  are  shown  in  Figure  6.14  and 
Table  6.  Since  the  experimental  data  points  for  the  base  pressure 
coefficients  were  not  obtained  at  the  same  pressure  ratio  values  as 
the  schlieren  data,  experimental  values  for  the  base  pressure  coeffi¬ 
cients  in  Table  6  were  interpolated  from  the  available  data  points 
at  the  five  given  nozzle  pressure  ratios.  Values  of  the  computed 
nozzle  base  pressure  were  in  good  agreement  with  the  experimental 
data  (3-7 7  error),  with  the  exception  of  the  highly  overexpanded 
case  at  which  P^/P^  =  0.150.  Figure  38  shows  that  as  the  pressure 
ratio  of  the  nozzle  is  lowered,  the  decreasing  trend  in  nozzle  base 
pressure  reverses  at  a  value  of  approximately  P^/P^  “  0.18  and  sharply 
increases  as  the  pressure  ratio  is  further  reduced.  This  sudden  reversal 
in  behavior  is  apparently  due  to  flow  separation  in  the  divergent 
portion  of  the  nozzle  which  prevents  the  jet  flow  from  expanding 
fully  to  its  assumed  Mach  3.0  state.  For  pressure  ratio  values  less 
than  P./P  »  0.135,  deterioration  of  the  Mach  disc  formation  to  a 

j  oo 

regular  shock  reflection  occurs  as  shown  in  Figure  39.  This  may 
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O  KXPKR IMKNTAL  DATA  (16) 
-O-  COMPl'TATIONAI.  SOLITIONS 


P|/Pco  =  0.124  P]/Pcd  =  0.132 


Pj /Pod  =  0.135  Pj/ Pod  =0.150 


Figure  39.  Schlicren  Photographs  (16)  Showing  the  Eventual 
Deterioration  of  the  Mach  Disc  with  Decreasing 
Nozzle  Pressure  Ratio. 
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indicate  that  for  values  less  than  P . /P  =  0.135,  either  the  jet 
Mach  number  is  less  than  1.48  and  thus  cannot  support  a  Mach  disc 
structure  (8) ,  or  the  increase  in  nozzle  static  pressure  due  to  the 
reduced  expansion  cannot  produce  the  deflection  angle  in  the  jet  flow 
needed  for  a  Mach  disc  to  occur.  A  non-separated  condition  was 
assumed  by  the  experimental  investigators,  since  their  value  of  P 
was  determined  using  the  jet  total  pressure  and  the  final  area  ratio 
of  the  nozzle.  Likewise,  the  computational  solutions  assumed  non- 
separated  Mach  3  flow  just  upstream  of  the  nozzle  exit  plane.  If  some 
separation  did  occur  and  the  nozzle  flow  did  not  fully  expand  to  a 
Mach  3  condition,  a  substantial  difference  in  base  pressure  could 
result . 

This  hypothesis  of  separation  in  the  nozzle  was  partially  confirm¬ 
ed  by  computationally  solving  a  case  where  the  jet  total  pressure 
corresponded  to  the  attached  case  (P^/P^  =  0.150),  but  with  jet 
input  profiles  corresponding  to  an  isentropic  expansion  cf  the  jet 
flow  to  a  Mach  number  of  only  2.60  as  detailed  in  Appendix  E.  A 
correct  strong  shock  structure  was  obtained  computationally,  and 
the  value  of  the  base  pressure  coefficient  increased  to  a  value  of 
-0.265.  This  was  in  much  better  agreement  with  the  experimental 
data  at  this  condition.  A  more  accurate  simulation  of  this  pressure 
ratio  condition  would  require  extending  the  computational  mesh  back 
to  the  nozzle  throat.  Separation  in  the  nozzle  could  then  occur  in 
a  direct  manner  in  the  numerical  solution.  Since  this  would  require 
extensive  grid  revisions  as  well  as  additional  computer  resources, 
it  was  considered  to  be  beyond  the  scope  of  this  investigation. 
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CHAPTER  VII 


CONCLUSIONS  AND  RECOMMENDATIONS 

A  numerical  method  of  obtaining  solutions  to  the  Navier-Stokes 
equations  for  supersonic  coflowing  axisymmetric  nozzles  has  been 
developed  from  a  selection  of  appropriate  techniques.  Based  on  the 
numerical  analysis  and  computational  results  obtained  through  this 
study,  the  following  conclusions  are  drawn: 

1.  The  numerical  solution  of  the  Navier-Stokes  equations  applied  to 
supersonic  coflowing  nozzles  successfully  reproduced  all  of  the 
essential  flow  features  including  boundary  layers,  corner  expan¬ 
sions,  recompression  shocks,  the  recirculation  region  adjacent 
to  the  nozzle  base  wall  and  the  evolution  of  the  near  wake  to  a 
flow  with  far  wake  behavior. 

2.  The  numerical  method  achieved  a  correct  transition  from  regularly 
reflected  shock  waves  at  the  line  of  symmetry  in  the  jet  core 
flow  to  a  strong  Mach  disc  formation  at  the  appropriate  static 
pressure  ratio  condition  of  the  nozzle.  The  subsonic  embedded 
region  immediately  behind  the  Mach  disc  formation  was  simulated 
in  a  correct  manner. 

3.  The  application  of  an  adaptive  grid  scheme  in  the  wake  region 
of  the  nozzle  annulus  successfully  positioned  the  fine  mesh 
region  of  the  computational  grid  in  the  wake  region  which  nor¬ 
mally  contains  severe  flow  gradients.  This  allowed  the  accurate 
simulation  of  this  high  flow  gradient  region  while  conserving 


numerical  resources. 


4.  The  nozzle  base  pressure  was  heavily  dependent  on  the  eddy 

viscosity  model  applied  in  the  region  of  the  near  wake.  Once 

the  model  was  tuned  for  the  neutrally  expanded  case  (P./P  = 

J  “ 

1.03),  good  agreement  was  obtained  computationally  for  all  cases 
where  the  flow  obeyed  the  assumption  of  remaining  attached  in 
the  divergent  portion  of  the  nozzle. 

5.  Boundary  conditions  must  be  carefully  formulated  and  applied 

in  order  to  prevent  physically  unrealistic  results  or  numerical 
divergence  of  the  solution.  Both  the  centerline  and  downstream 
boundaries  were  sensitive  where  regions  of  subsonic  flow  occured 
over  a  substantial  portion  of  the  boundary.  Both  the  quadratic 
extrapolation  used  in  regions  of  supersonic  flow  and  a  second  order, 
zero  gradient  condition  caused  either  divergence  or  unrealistic 
conditions  at  the  boundary  when  applied  to  regions  of  subsonic 
flow.  A  first  order  zero  gradient  condition  was  used  successfully 
in  these  regions  of  subsonic  flow  and  found  to  be  superior. 
Generalizations  about  the  success  of  this  first  order  boundary 
condition  cannot  be  made,  since  the  degree  of  success  achieved 
is  dependent  on  the  specific  numerical  algorithm  applied  to  the 
problem. 

6.  The  final  steady  state  solutions  were  found  to  be  insensitive  to 
the  initial  conditions  applied  over  the  computational  domain. 
However,  the  time  to  converge  to  the  final  solution  was  highly 
dependent  on  the  application  of  specific  initial  conditions.  In 
particular,  the  region  of  subsonic  recirculation  in  the  near 
wake  was  the  last  region  in  the  solution  domain  to  converge. 
Solutions  started  using  only  the  boundary  layer  profiles  across 
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the  domain  required  three  hours  to  converge  on  a  Cyber  175  (45  x 
45  point  mesh),  whereas  cases  started  from  a  previous  different 
jet  total  pressure  condition  but  with  an  established  near  wake 
structure  required  only  1.7  hours  to  converge  (45  x  45  mesh). 

7.  To  the  author's  knowledge,  this  is  the  first  full  Navier-Stokes 

solution  that  has  accurately  simulated  the  viscous-inviscid  inter¬ 
actions  present  in  a  supersonic  coflowing  nozzle  at  off-design 
conditions  where  the  strong  Mach  disc  shock  structure  is  present. 
Mikhail  (8)  previously  was  unsuccessful  in  reproducing  the  Mach 
disc  reflection  in  a  full  Navier-Stokes  solution  due  to  the  prob- 
bable  improper  placement  of  the  jet  inflow  boundary  condition, 
which  did  not  allow  the  jet  plume  to  expand  to  the  degree  necessary 
to  generate  a  Mach  disc  reflection. 

Based  on  the  numerical  analysis  and  results  obtained  through 

this  study,  the  following  recommendations  are  made: 

1.  The  present  scalar  computer  code  developed  during  the  course  of 
this  investigation  should  be  vectorized  for  use  on  the  new  genera¬ 
tion  of  "supercomputers"  such  as  the  CRAY-1  or  the  Cyber  203. 
Although  present  solution  times  are  on  the  order  of  two  to  four 
hours  when  run  on  a  Cyber  175  computer,  a  fully  vectorized  version 
of  the  present  computer  code  can  be  expected  to  converge  within 
five  minutes  on  a  CRAY-1  computer  (39).  This  will  allow  computa¬ 
tion  of  more  complex  nozzle  geometries  and  better  resolution  in 
the  boundary  layers  through  the  application  of  finer  mesh,  while 
holding  costs  to  a  reasonable  level. 

2.  The  present  numerical  solver  should  be  modified  to  include  the 
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effects  of  a  calorically  imperfect  gas  with  variable  specific  heat 
and  thermal  conductivity.  This  modification  would  allow  accurate 


simulation  of  the  temperature  dependent  effects  for  hot  exhaust 
nozzles  with  gas  temperatures  less  than  5000°R  (19).  Only  minor 
revisions  to  the  existing  computer  code  would  be  required  in  order 
to  include  these  temperature  effects. 

3.  After  the  implementation  of  the  two  previous  recommendations, 

it  would  be  desirable  to  incorporate  the  effects  of  species  mixing 
into  the  numerical  solver.  Many  practical  cases  of  interest  involve 
a  jet  exhaust  gas  with  a  different  species  than  that  of  the  external 
stream.  This  modification  would  require  a  significant  code  revision, 
since  the  addition  of  the  equation  of  mass  diffusion  would  be  re¬ 
quired,  as  would  the  correct  modeling  of  appropriate  mass  diffusion 


coef  f ic ients  . 
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APPENDIX  A 

NOZZLE  WALL  TEMPERATURE  CALCULATION 


The  nozzle  wall  boundary  condition  applied  during  the  numerical 
solution  procedure  assumes  a  constant  nozzle  wall  temperature  within 
the  computational  domain.  The  relatively  high  thermal  conductivity  of 
the  stainless  steel  nozzle  makes  this  assumption  valid.  This  wall 
temperature  can  be  calculated  by  applying  a  heat  flux  balance  across 
the  freestream  boundary  layer,  the  nozzle  wall,  and  the  jet  boundary 
layer  as  shown  in  Figure  40.  Conduction  of  heat  in  the  axial  direction 
is  neglected  due  to  the  low  temperature  gradients  in  this  direction. 

Since  both  the  freestream  and  the  jet  flow  are  of  a  high-velocity 
nature,  the  unit  heat  flux  for  either  stream  can  be  written  as  (40): 


q  =  h.  (T  -  T  ) 
nw .  1  aw .  w 

l  l 


(A-l) 


where  h.  is  the  heat  conductance  of  each  flowstream.  The  adiabatic 

l 

wall  temperature,  T  ,  is  defined  by  the  expression: 

3W  . 

1 


T  =  T  (1  +  n/p7  M  2) 

aw ,x  2  i 

l 


(A-2) 

where  n  is  1/2  for  laminar  flow  and  1/3  for  turbulent  flow.  Equating 

the  heat  flow  out  of  the  control  volume  for  a  steady  state  process  gives 

z  q,  -  0  (A-3) 

i 


2n  rJ„°  q  +  2n  r^°  q.  =  0 
J  J 


(A-4) 


where  j0  is  either  0  or  1  for  two-dimensional  or  axisymmetric  flow, 
respectively . 
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Figure  40.  Heat  Flux  Balance  Used  to  Determine  the 
Nozzle  Wall  Temperature. 
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Applying  equation  A-l  to  the  values  of  the  individual  heat  fluxes  and 

solving  for  the  wall  temperature  gives: 

T  +  (r0/r.)j°(h  /h . ) T 
aw .  1  00  j  aw 

T  = - J - — 7 — -  (A- 5) 

W  [1  +  <ra/r  J^tti^/h  )1 

For  flows  involving  a  constant  boundary  layer  edge  velocity  and  a 
constant  wall  temperature,  the  conductances,  h^,  can  be  defined  by  the 
following  analytic  expressions: 

-2/3 

h.  =  0.332  Cn  p.  u.  Re  ~  Pr  '  (A-6) 

1  P  1  1  xi 

for  laminar  flow,  and 

h.  =  0.0295  Cn  p.  u.  RS1^5  Pr_2/5  (A-7) 

1  P  1  X  XI 

for  turbulent  flow  (40).  All  of  the  quantities  on  the  right-hand-side 
of  equation  A-5  are  then  known,  and  a  value  for  the  wall  temperature 
can  be  calculated  based  on  the  states  of  the  two  flowfields  adjacent 
to  the  nozzle  wall. 
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APPENDIX  B 

ADAPTIVE  FINITE  DIFFERENCE  MESH 


It  is  desirable  that  the  fine  mesh  region  of  the  computational 
grid  remain  in  the  areas  of  relatively  high  velocity  and  temperature 
gradients  as  the  solution  progresses  towards  convergence.  Hirt  (Tl) 
has  used  a  technique  in  the  solution  of  free  surface  flows  that  allows 
the  grid  to  adapt  as  the  solution  progresses.  The  following  kinematic 
equation  is  applied  in  the  region  where  the  nozzle  wake  and  shear  layer 
develop : 

f  =  CA(v  -  u  f  >  <B-1) 

This  equation  ensures  the  condition  that  as  the  solution  converges, 
the  physical  slope  of  the  constant  n  finite  difference  cell  boundaries 
is  the  same  as  that  of  the  velocity  vectors  near  each  cell. 

Equation  ( B— 1 )  can  be  converted  to  the  following  finite  difference 
form  for  application  to  a  computational  mesh: 


rn+1  =  rn  .  +  CA  At  [  v?  .  -  u" 

1,1  1,1  A  i,j  i,j  n  n 

x.  .  -  x.  ,  . 

1,1  t-1, J 


n  n 

r .  .  -  r .  .  . 
i,l  i-l,  h 


where 


n 

+  un 

+  un 

+  un 

i-l.J 

i+1  ,j 

i  ,J-1 

i  ,  j  +  1 

n 

+  vn 

+  vn 

+  vn 

i-l>j 

i+1 ,  j 

i ,  j-1 

i,J  +  l 

and  where  is  a  constant  used  to  damp  the  grid  motion  with  respect 
to  time.  Spatial  averaging  of  the  velocity  components  is  applied  in 
order  to  reduce  the  effects  of  numerical  velocity  fluctuations  at  indi- 
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vidual  mesh  points  in  the  flow.  The  upwind  difference  form  of  the 
cell  aspect  ratio  term  (Ar/Ax)  is  also  utilized  to  achieve  better 
stability  in  the  finite  difference  equation. 

Equation  (B-l)  is  applied  in  the  wake  region  of  the  flowfield 
for  a  line  of  constant  n  (j=constant) ,  where  the  specific  value  of  n 
corresponds  to  the  nozzle  inner  wall  for  overexpanded  flowfields,  and 
to  the  nozzle  outer  wall  for  underexpanded  flowfields.  Once  the  posi¬ 
tion  of  this  grid  reference  line  is  established,  the  fine  mesh  region 
corresponding  to  the  wall  thickness  in  computed.  The  exponential 
stretching  scheme  discussed  in  Chapter  III  is  then  applied  for  each 
value  of  J,  in  the  regions  above  and  below  the  fine  mesh  region  as  shown 
in  Figure  41.  The  first  two  grid  points  above  the  centerline  were  kept 
fixed  at  constant  heights  for  all  values  of  x.  This  prevented  large 
numerical  errors  in  axisymmetric  cases  involving  the  differencing  of 
terms  containing  (1/r),  where  r  is  a  very  small  number. 

The  constant  was  specified  in  the  range  of  0.3  -  0.6  in 
order  to  allow  the  grid  to  adapt  smoothly  as  the  solution  converged. 
Larger  values  of  caused  undesirable  oscillatory  motion  of  the  grid 
reference  line  with  respect  to  time. 

The  adaptive  grid  scheme  was  applied  once  during  every  iteration 
of  the  solution  algorithm.  The  number  of  points  allowed  to  "float"  on 
the  grid  reference  line  could  be  varied  during  the  course  of  the  solu¬ 
tion.  This  capability  was  utilized  primarily  during  the  start-up 
portion  of  the  numerical  solution,  where  only  a  limited  number  of  points 
close  to  the  nozzle  were  allowed  to  "float"  until  the  shear  layer  was 
established.  After  the  position  of  the  shear  layer  across  the  complete 
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computational  domain  was  completely  established,  the  adaptive  grid 
scheme  could  be  turned  off  in  order  to  save  computer  time  during  the 
remainder  of  the  solution. 


113 


APPENDIX  C 


TWO-DIMENSIONAL  FLAT  PLATE  FAR  WAKE  SOLDI  ION 

In  order  to  test  the  validity  of  the  turbulence  model  in  the 
far  wake  region,  a  check  case  possessing  known  experimental  data  in 
this  region  was  computed.  The  exterimental  data  for  this  case  is 
that  of  Toyoda  and  Hiriyama  (33),  which  involved  a  two-dimensional 
thin  flat  plate  at  a  Mach  number  of  1.60.  The  flat  plate  possessed 
a  thickness  of  1  mm  and  a  trailing  edge  thickness  of  0.1  mm. 

The  following  flow  quantities  were  given  in  the  experimental 

data : 

M  =1.60 

oo 

P0  =  3.2  atm 

<X 

Re  =  5300 

A  value  for  free  stream  stagnation  temperature  was  not  given  experi¬ 
mentally.  A  stagnation  temperature  of  518. 7°R  was  therefore 
assumed.  Lsing  the  previous  quantities,  the  value  for  the  boundary 
layer  momentum  thickness  at  the  trailing  edge  was  computed  as: 

=  . 0895  mm 

The  computational  solution  used  the  same  numerical  silver 
(Macf.ormack's  explicit  method)  as  that  of  the  coflowing  nozzle  s,  lut 
The  eddy  viscosity  models  utilized  were  identical  with  those  used  in 
the  boundary  layer  and  far  wake  regions  of  the  coflowing  nozzle .  ! 

boundary  conditions  utilized  also  closely  resemble  thost  ot  tit  .of! 
nozzle.  The  inflow  conditions  were  set  by  the  experimental  data  and 
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remained  fixed  for  the  duration  of  the  solution.  The  upstream  u 


component 

of  velocity  near  the  trailing  edge  of  the 

plate 

was  matched 

to  that  of 

d it  ions  at 

the  experimental  data  as  shown  in  Figure 

the  upstream  boundary  were  then  set  as: 

42. 

Input  con- 

u  ( x  .  , 

l 

y,  t)  =  0.6295  Uoo(y/uTF)1/5 

(C-l) 

v(x.  , 

y,  t)  =  0.0 

(C-2) 

‘  (v 

y,  t)  =  .rx/{l+  ^  np— ^M“  [  1  -  (u(x.,  y,  t )  / 

0 

uj‘1 

}  (C—  3) 

P(x.  , 

y ,  t )  =  Pi 

(C-4) 

Both  the  upper  and  lower  freestream  boundaries  utilize  the  charact¬ 
eristic  condition  applied  at  the  upper  boundary  of  the  coflowing 
nozzle  solutions.  The  outflow  and  wall  boundaries  used  conditions 
identical  to  those  utilized  for  the  coflowing  nozzle.  The  computational 
solution  was  initially  started  by  applying  the  upstream  profile  across 
the  (orputation.il  domain: 

"  i  >: .  .  t*>  =  "  i  >:  j  ,  y,  o)  (C— 5) 

1 1  i  1  a  I .  ulat  ion  ut  ilized  a  3^x34  computational  mesh,  with  ex- 
f  •  ■  •  ■  -  - :  t  ;  a  1  r.e-d  it  reti  King  employed  in  both  the  x  and  y  directions 
as  s'  own  in  Henri  43.  Minimum  grid  spacings  in  the  x  and  y  direc¬ 
ti  n-  wert  Hi  2  Jr:-.  and  0 . 1  OOnm ,  respectively.  The  physical  dimensions 
of  the  cor.putu: ional  flowfield  n  the  x  and  y  directions  were  38mra 
and  Jinn.,  respectively.  The  rati  of  data  processing  on  a  CDC  Cvber 
173  computer  wan  f'.0O14  sec  per  grid  point  for  each  iterative  time 
step.  The  solution  was  computed  for  a  duration  of  four  characteristic 
times  ( 3d 00  iterations),  at  which  tine  no  significant  change  was 
detected  in  the  dependent  variables.  The  result  was  then  taken  to 


y  /e 


lit’urr  A'l.  Velocity  Profile  at  the  Trailing  I.dpe  of 
the  Two-Dimensional  Flat  Plate. 
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be  the  asymptotic  solution. 


The  computational  solution  demonstrated  that  the  eddy  viscosity 
model  gives  good  results  in  the  far  wake  region.  The  maximum  velocity 
defect  generated  computationally  is  in  good  agreement  with  the  ex¬ 
perimental  data  as  shown  in  Figure  21.  Figure  43  shows  that  the 
velocity  fic1d  generated  computationally  evolves  from  the  boundary 
layers  on  the  plate  to  a  classic  wake  solution  very  rapidly  due  to 
the  turbulent  nature  of  the  flow. 
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APPENDIX  D 


TWO-DIMENSIONAL  WEDGE-FLAT  PLATE  NEAR  WAKE  SOLUTION 

In  order  to  test  the  validity  of  the  turbulence  model  utilized 
in  the  near  wake  region  of  the  coflowing  nozzle  solutions,  a  check 
case  exhibiting  similar  physical  characteristics  and  possessing 
known  experimental  data  in  this  region  was  solved  numerically.  The 
experimental  case  selected  for  this  validation  study  is  one  obtained 
by  Rom,  Seginer  and  Kronzon  (35)  for  a  two-dimensional  wedge-flat 
plate  in  a  turbulent  supersonic  flowfield.  The  model  used  for  this 
study  consisted  of  a  sharp  15°  half  angle  wedge-flat  plate  with  a 
base  height  of  iO  mm  and  a  chord  of  44  mm.  The  following  flow  con¬ 
ditions  were  given  in  the  experimental  data: 

M  =2.25 

03 

Re  =  1.5  x  106 
c 

P  =40  psig 

00 

T  =  492°  R 
o 

CO 

A  computed  adiabatic  wall  temperature  of  466.6°  R  based  on  the  flow 
condition  at  the  flat  plate  portion  of  the  model  was  used  in  the 
computational  solution. 

The  numerical  solution  used  the  same  computational  solver  and 
turbulence  model  as  that  of  the  coflowing  nozzle  solutions.  The 
boundary  and  initial  conditions  were  also  identical  to  those  used 
in  the  coflowing  nozzle  with  the  exception  of  the  jet  centerline 
condition.  As  shown  in  Figure  45,  this  condition  was  replaced  with 
a  lower  freestrcam  boundary  which  utilized  a  characteristic  scheme 
similar  to  that  of  the  upper  freestream  boundary.  Since  a  value  for 
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gure  45.  Computational  Boundary  Conditions  for  the  Two-Dimensional 
Wedge-Flat  Plate. 


the  boundary  layer  thickness  near  the  trailing  edge  was  given  experi¬ 
mentally  (i/h  =  0.15),  a  two-dimensional  boundary  layer  code  was 
used  to  generate  the  velocity  and  temperature  profiles  on  the  up¬ 
stream  boundary.  The  boundary  layer  starting  length  was  adjusted  to 
give  the  correct  boundary  layer  thickness.  One  additional  condition 
was  imposed  on  the  line  of  symmetry  for  the  wedge  (j  =  23) .  At  this 
line  of  symmetry,  a  zero  v  velocity  component  was  enforced  to  help 
stabilize  the  wake  during  the  startup  of  the  numerical  solution  and 
help  accelerate  convergence  by  damping  any  numerical  shear  layer 
oscillations  in  the  wake. 

The  solution  was  calculated  using  a  45  x  45  point  computational 
mesh  with  exponential  stretching  employed  in  both  the  x  and  y 
directions  as  shown  in  Figure  46.  The  physical  dimensions  of  the 
computational  flowfield  in  the  x  and  y  directions  were  10cm  and  7cm, 
respectively.  Minimum  grid  spacing  in  the  x  and  y  directions  was 
set  at  0.5  mm.  This  gave  a  value  of  u/u^  =  0.83  for  the  first  point 
in  the  boundary  layer  above  the  nozzle  wall,  which  corresponds  to  an 
identical  value  in  the  jet  boundary  layer  of  the  coflowing  nozzle 
solutions.  Thus,  truncation  error  should  be  similar  in  this  region 
for  both  the  wedge-flat  plate  and  the  coflowing  nozzle.  The  rate 
of  data  processing  on  a  CDC  Cyber  175  computer  was  0.0014  sec  per 
grid  point  for  each  iterative  time  step.  The  solution  was  allowed 
to  progress  for  approximately  four  characteristic  timesteps  at  which 
time  the  change  in  the  dependent  variables  was  less  than  0.5%  per 
characteristic  time  period.  This  condition  was  then  considered  to 
be  the  converged  asymptotic  solution. 
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Specific  features  o.  the  experimental  flowfield  in  the  near 
wake  were  reproduced  in  the  computational  solution  and  can  be  dis¬ 
tinguished  in  the  plots  of  Mach  number  contours  and  velocity  profiles 
shown  in  Figures  47  and  48.  Several  flow  features  which  were  numeri¬ 
cally  observed  include  the  existence  of  the  boundary  layers  along 
thn  horizontal  walls  of  the  body,  the  turning  of  the  flow  through  the 
corner  expansion  fans,  the  existence  of  the  subsonic  recirculating 
"dead  water"  region  adjacent  to  the  base  of  the  body,  flow  recom¬ 
pression  through  the  trailing  shocks,  and  the  evolution  of  the  wake 
to  a  classic  far  wake  flow.  The  weak  lip  shock  evident  in  the 
experimental  data  was  not  readily  evident  in  the  numerical  solution. 
This  may  be  attricited  to  the  fact  that  the  numerical  method  tends 
to  smear  shocks,  and  thus  has  difficulty  locating  shocks  which  are 
very  weak. 

Quantitative  accuracy  of  the  numerical  solution  is  identified 
through  the  use  of  the  given  experimental  static  pressure  and  pitot 
pressure  data.  A  comparison  of  the  axial  static  pressure  distribu¬ 
tion  along  the  line  of  symmetry  is  shown  in  Figure  49.  The  computed 
static  pressure  distribution  is  within  57  of  the  measured  values 
except  in  the  region  of  recirculating  flow  (x  1.0cm),  where  there 
is  up  to  a  7%  discrepancy  between  the  experimental  and  computed  values. 
However,  static  pressure  probes  like  the  one  used  to  obtain  the  data 
are  very  sensitive  to  flow  angularity,  and  thus  less  leliable  in  areas 
of  recirculating  flow.  A  more  reliable  comparison  is  that  of  the 
pitot-pressure  surveys  shown  in  Figure  SO.  These  measurements  are 
less  susceptible  to  errors  resulting  from  flow  angularity.  This 
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igure  47.  Two-Dimensional  Wedge-Flat  Plate  With  the  Computed  Mach  Number 
Contours  Shown. 


Figure  48.  Computed  Velocity  Profiles  in  the  Near  Wake  of  the 
♦  Two-Dimensional  Wedge-Flat  Plate. 


i 

i  Figure  49.  Static  Pressure  Along  the  Line  of  Symmetry  in  the 

;  Near  Wake  of  the  Two-Dimensional  Wedge-Flat  Plate. 

I 

t 


^  1 28 


l 


figure  shows  excellent  agreement  in  the  "deadwa ter ”  region  of 
recirculating  flow.  The  numerical  solution  smears  the  weak  be¬ 
ginning  of  the  recompression  trailing  shocks  at  0.5cm,  but  correctly 
simulates  them  at  the  correct  values  further  downstream  (x  >  1.0). 
This  figure  particularly  demonstrates  that  the  phenomena  present  in 
the  near  wake  are  accurately  simulated  by  the  present  computational 
method  and  turbulence  modeling. 


t 
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APPENDIX  E 


AX I SYM' \ ET P. T C.  CO F LOWING  NO Z Z_L E  SOLUTION 
S IMPUTING  INTERNAL  SEPARATION  OF_THE  NOZZLE 

The  peculiar  reversal  in  the  nozzle  base  pressure  coefficient 

as  the  pressure  ratio  was  lowered  below  P./P  =  0.18  was  bilieved  to 

J  °" 

be  caused  by  separation  in  the  nozzle  (Fig.  38).  This  was  further  sub¬ 
stantiated  through  an  examination  of  the  schlieren  photographs  (16) 
which  show  a  definite  change  in  flow  character  near  this  value  (Fig.  39). 
To  test  this  hypothesis  a  numerical  solution  was  computed  for  a  case 
where  the  jet  total  pressure  corresponded  to  that  of  the  attached  strong 

shock  solution  (P./P  =  0.15),  but  where  the  jet  Mach  number  equals 

J  00 

2.60  based  upon  an  isentropic  expansion  in  the  nozzle.  For  this  jet 
Mach  number  a  reduction  in  nozzle  area  ratio  was  assumed  from  A/A*  =  4.2 
at  Mach  3  to  A/A*  =  2.9  at  Mach  2.60.  This  reduced  expansion  rate 
was  intended  to  roughly  simulate  the  effects  of  boundary  layer  separa¬ 
tion  in  the  nozzle,  while  retaining  grid  geometry  and  fineness  identical 
to  the  previously  calculated  attached  jet  flow  case. 

The  following  jet  flow  parameters  were  used  in  this  solution: 

M.  =  2.60 
J 

T  .  =  580. 5°R 
°J 

P  .  =  1636  psf 

Oj 

These  jet  conditions  produce  an  actual  nozzle  pressure  ratio  of 

P  /P  =  0.276  versus  the  calculated  value  of  P./P_  =  0.150  which  was 
J  "•  J 

assumed  in  P.eference  16. 

This  solution  was  initialized  using  the  solution  for  the  attached 
strong  shock  case  with  a  calculated  nozzle  pressure  ratio  equal  to 
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.  i  '•  .  !  ■:>  »■;  t  !ur  the  upst  rear;  jet  boundary,  boundary  conditions 

.t:.d  -  i  '.  i  i'll!  i^urat  ion  identical  to  those  of  the  attached  case  were 
utilised.  At  the  upstream  jet  boundary  the  following  velocity  pro- 


!  ili  wa  :  i  xed  : 


0  ■  r/r.  -  0.835 
jet  - 

0.83:  <  r/r.  <  1 
jet  - 

0  <  r/r.  _  •  1 
jet  - 


u  =  2002  ft/sec 

u  =  0 

v  =  0 

As  shown  in  Figure  51,  this  velocity  input  profile  increased 

the  displacement  thickness  of  the  jet  boundary  layer  in  rough  approxi¬ 
mation  to  that  caused  by  separation  in  the  nozzle.  In  this  case  the 
jet  flow  is  assumed  to  be  turbulent  in  nature. 

As  shown  in  Figure  52,  a  strong  shock  solution  was  obtained,  with 

a  reflection  length  of  x  lx.  -  1.29,  about  2%  greater  than  the  ex- 
e  s  jet 

perimental  value  at  this  jet  total  pressure  case.  The  base  pressure 
coefficient  obtained  was  equal  to  -0.264,  about  8Z  less  than  the 
experimental  value  and  in  much  better  agreement  than  that  of  the 
attached  jet  boundary  layer  case. 

Although  discrepancies  in  this  solution  such  as  the  smaller 
diameter  of  the  computational  Mach  disc  and  the  appearance  of  a 
physically  nonexistant  secondary  Mach  disc  further  donwstrean  are 
apparent,  this  case  does  demonstrate  that  separation  and  its  effect 
on  the  flow  expansion  in  the  nozzle  can  significantly  impact  the 
resultant  base  pressure  coefficient  of  the  nozzle.  Therefore, 
future  computations  for  low  pressure  ratios  should  commence  at  the 
throat  section  to  insure  that  nozzle  separation  is  properly  considered. 
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the  Separated  Flow  Simulation. 


APPENDIX  F 


1NVESTIGAT ION  OF  N UMER I  CAL  ERROR 

Two  error  analyses  were  performed  on  the  numerical  technique 
used  to  obtain  solutions  for  the  coflowing  nozzle.  The  first  con¬ 
sisted  of  examining  the  error  generated  by  the  numerical  algorithm 
used  to  solve  the  Navier-Stok.es  equations.  The  second  was  a  study 
of  the  effect  of  repositioning  first  the  downstream  boundary,  and 
subsequently  the  upper  boundary  to  regions  containing  only  minor  flow 
gradients  normal  to  each  boundary.  These  analyses  are  discussed  in 
detail  in  the  following  sections. 


Truncation  Error  Analysis 

MacCormack's  explicit  finite  difference  algorithm  is  an  equiva¬ 
lent  second  order  accurate  numerical  solver.  The  final  converged 
solution  for  any  case  computed  by  this  algorithm  should  then  satisfy 
the  Navier-Stokes  equations  at  all  interior  node  points  with  second 
order  accuracy.  A  numerical  check  on  this  accuracy  was  conducted 
using  a  typical  converged  solution  and  the  following  procedure. 

As  shown  in  Chapter  2,  the  axisymmetric  Navier-Stokes  equations 
can  he  written  as: 


:  +  +  -1  -  »  =  o 

it  >x  r  Jr  r 


(F-l) 


A  nond imens i onal ized  finite  difference  formulation  of  these  equations 
can  then  be  written  as: 


(-c)  *1  +  £  [M  +  I  Aii£>  _ 

D  At  V  Ax  r  Ar  rJ 


E  =  Error 


(F-2) 
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where  t’  is  defined  as: 


(F— 3) 


and  t  =  5  x  10  sec  is  set  to  insure  that  the  lead  terms  on  the 
c 

Teft  hand  side  of  equation  ( F— 2 )  are  of  order  one. 

The  particular  check  case  selected  was  that  containing  the 
strong  Mach  disc  shock  structure  (P^/P^  =  0.150),  since  this 
case  contains  both  substantial  regions  of  subsonic  and  supersonic 
flow,  and  oblique  as  well  as  normal  shock  waves.  The  MacCormack 
solution  for  this  case  was  used  as  input  for  equation  (F-3),  where 
the  left  hand  side  of  the  equation  was  computed  at  all  interior 
grid  points  using  a  standard  two-dimensional  second  order  central 
differencing  scheme  applied  on  the  transformed  computational  plane. 
The  magnitude  of  the  Error  vector  (E)  is  then  an  indication  of  how 
close  MacTornack's  method  is  to  an  alternate  second  order  accurate 
solver.  The  following  root  mean  square  (RMS)  values  of  E  were 
obtained  over  the  interior  of  the  computational  domain: 


f  E 

p 

> 

0.026 

E 

rms 

'  f: 

pu 

= 

0.020 

E 

PV 

0.011 

E 

l  pe 

0.026 
<  / 

This  result  indicates  that  over  the  domain,  MacCormack's  algorithm 
and  the  two-dimension;  .di  difference  scheme  are  equivalent  to 

within  three  percent. 
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Finally,  the  Error  vector  E  was  examined  over  the  computational 
domain  to  determine  which  regions  generated  the  highest  magnitude  of 
error.  Since  the  RMS  error  value  of  the  continuity  equation  was  one 
of  the  largest,  and  the  error  distribution  was  representative  of  that 
in  the  other  equations,  it  is  shovm  in  Figure  53.  In  this  figure, 
only  error  values  greater  than  the  RMS  value  are  shown  as  contours. 
Regions  containing  the  largest  error  consist  of  those  containing  shock 
waves  and  that  containing  the  expansion  fan  near  the  sharp  corner 
of  the  nozzle.  In  these  shock  regions,  strong  flow  gradients  exist 
over  areas  with  fairly  coarse  finite  difference  mesh  spacing.  Although 
the  wake  region  also  contains  strong  gradients  within  the  mixing  layer, 
it  lies  within  a  fine  mesh  region  of  the  grid  that  produces  much  less 
numerical  error. 

This  analysis  further  demonstrates  the  desirability  of  utilizing 
adaptive  mesh  schemes  that  can  align  the  grid  with  flow  gradients 
as  the  solution  progresses  to  convergence. 

Boundary  Position  Analysis 

Although  the  upstream  boundary,  the  centerline  boundary,  and 
the  position  of  the  nozzle  walls  were  fixed  by  the  definition  of  the 
problem  to  be  solved,  the  placement  of  both  the  downstream  boundary 
and  the  upper  boundary  was  left  to  the  discretion  of  the  computational 
investigator.  It  was  desirable  to  place  these  boundaries  as  close  to 
the  nozzle  as  possible  in  order  to  achieve  computational  efficiency 
in  a  compact  domain.  In  the  coflowing  nozzle  solutions,  both  of  these 
boundaries  were  located  in  regions  in  which  flow  gradients  existed  due 
to  the  presence  of  shock  and  expansion  waves  and  viscous  phenomena  such 
as  shear  layers  and  wakes  in  the  flowfield.  The  assumption  is  made  that 
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MACH  NUMBER  CONTOURS 


Figure  53.  Comparison  Between  the  Computational  Solution 
and  the  Error  Present  in  the  Continuity 
Equation,  P./P  =  0.150. 


positioning  these  two  boundaries  in  flow  gradient  regions  does  not 
affect  the  computational  solutions  obtained.  To  validate  this 
assumption,  each  of  these  boundaries  was  repositioned  a  greater  dis¬ 
tance  from  the  nozzle  to  regions  containing  only  minor  flow  gradients 
normal  to  each  boundary.  The  resulting  effect  on  the  shock  wave 
structure  as  well  as  on  the  nozzle  base  pressure  coefficient  in  the 
numerical  solution  was  then  observed. 

The  coflowing  nozzle  case  for  which  P./P  =  0.251  was  examined  for 

J  00 

this  particular  study.  The  downstream  boundary  contains  primarily  super¬ 
sonic  outflow  with  an  embedded  wake  region  of  subsonic  outflow.  The 


downstream  boundary  was  extended  from  its  original  position  at  x/r^^  =  3.0 
to  a  new  value  of  x/r^et  =  6.0.  This  stretching  was  achieved  by  the 
addition  of  twelve  grid  points  axially  to  the  original  mesh.  The  out¬ 


flow  at  this  new  position  was  totally  supersonic  in  nature  with  only 
minor  gradients  in  existence  normal  to  the  boundary.  As  shown  in 
Figure  54,  no  changes  were  evident  in  the  shock  structure  contained  in 


the  original  domain.  The  computational  base  pressure  coefficient 
remained  unchanged  at  a  value  of  p  =  -0.299. 

D 

The  effect  of  repositioning  the  upper  boundary  to  a  radial  distance 


at  which  flow  gradients  are  not  present  was  then  examined.  The  upper 

boundary  of  the  axially  stretched  case  was  extended  from  its  original 

value  of  r/r.  =  3.0  to  a  new  value  of  r/r.  =  6.0  as  shown  in  Figure 
jet  jet  0 

55.  In  this  case  eight  grid  points  were  added  radially  to  the  axially 
stretched  mesh.  Again,  no  changes  were  evident  in  the  numerical  shock 


structure,  and  the  computational  value  of  the  nozzle  base  pressure  co¬ 


efficient  remained  at  p  =  -0.299. 
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Since  repositioning  these  boundaries  had  essentially  no 
effect  on  the  computational  solution  that  was  examined,  the  applica¬ 
tion  of  these  boundary  conditions  in  the  original  regions  containing 
mixed  supersonic-subsonic  flow  on  the  outflow  boundary  and  substantial 
flow  gradients  on  both  boundaries  was  valid. 
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